diff --git a/ChangeLog b/ChangeLog index 39452ae74..82fa4118c 100644 --- a/ChangeLog +++ b/ChangeLog @@ -40,6 +40,28 @@ FlexibleSUSY 2.0.0 [not released yet] Thanks to Jobst Ziebell. + * Feature: BSM contributions to the electric dipole moment of + fermions at the 1L level in any given BSM model. Note: Diagrams + with non-SM vector bosons are not taken into account. + + In order to let FlexibleSUSY calculate the EDM of a fermion F, the + symbol FlexibleSUSYObservable`EDM[F] must be written into an SLHA + output block in the ExtraSLHAOutputBlocks variable in the + FlexibleSUSY model file. If F is a multiplet, the field index must + be specified, for example FlexibleSUSYObservable`EDM[F[1]] for the + first field in the multiplet. + + Example: + + ExtraSLHAOutputBlocks = { + {FlexibleSUSYLowEnergy, + {{23, FlexibleSUSYObservable`EDM[Fe[1]]}, + {24, FlexibleSUSYObservable`EDM[Fe[2]]}, + {25, FlexibleSUSYObservable`EDM[Fe[3]]} } } + }; + + Thanks to Jobst Ziebell. + * Feature: New functions, FSGetProblems[], FSGetWarnings[] and FSToSLHA[], have been added to FlexibleSUSY's spectrum generator Mathematica interface. The first diff --git a/doc/librarylink.dox b/doc/librarylink.dox index 433a2e569..6cb97baf2 100644 --- a/doc/librarylink.dox +++ b/doc/librarylink.dox @@ -227,7 +227,7 @@ _Output_: ~~~~~~~~~~~~~~~~~~~~{.m} { alphaEmMZ -> 0.00781763, (* alpha_em(MZ) in the SM(5), MS-bar *) - GF -> 0.0000116637, (* Fermi constant *) + GF -> 0.000011663787, (* Fermi constant *) alphaSMZ -> 0.1184, (* alpha_s(MZ) in the SM(5), MS-bar *) MZ -> 91.1876, (* Z pole mass *) mbmb -> 4.18, (* MS-bar bottom mass at Q = mb *) diff --git a/meta/CConversion.m b/meta/CConversion.m index 375e0c9bf..c11fc4630 100644 --- a/meta/CConversion.m +++ b/meta/CConversion.m @@ -237,7 +237,7 @@ CConversion`ArrayType[ CConversion`complexScalarCType,_] | CConversion`MatrixType[CConversion`complexScalarCType,__]| CConversion`TensorType[CConversion`complexScalarCType,__], - "(" <> expr <> ").cast<" <> CreateCType[CConversion`ScalarType[complexScalarCType]] <> " >()" + "(" <> expr <> ").template cast<" <> CreateCType[CConversion`ScalarType[complexScalarCType]] <> " >()" , _, Print["Error: CastTo: cannot cast expression ", expr, " to ", toType]; diff --git a/meta/EDM.m b/meta/EDM.m new file mode 100644 index 000000000..82641081c --- /dev/null +++ b/meta/EDM.m @@ -0,0 +1,540 @@ +(* ::Package:: *) + +BeginPackage["EDM`", {"SARAH`", "TextFormatting`", "TreeMasses`", "Vertices`", "Parameters`"}]; + +(* This module generates c++ code that calculates electric dipole moments of fields *) + +EDMInitialize::usage="EDMInitialize the EDM module."; + +EDMSetEDMFields::usage="Set the fields for which the EDMs shall be calculated."; +EDMCreateFields::usage="Returns the c++ code that contains all field fields"; +EDMCreateDiagrams::usage="Returns the c++ code that contains all relevant diagram classes"; +EDMCreateInterfaceFunctions::usage="Returns the c++ code containing the interface functions {prototypeCode, definitionCode}." +EDMCreateVertexFunctionData::usage="Returns the c++ code that contains all relevant vertex function data"; +EDMCreateDefinitions::usage="Returns the c++ that contains all function definitions" + +EDMNPointFunctions::usage="Returns a list of all n point functions that are needed. Actually it is a list of fake functions to extract vertex functions..."; + +(******** TODO: IMPORTANT NOTES: + If you add new kinds of vertices (e.g for new diagram types): + - Add the new types to vertexTypes + - Expand CouplingsForFields[] and VertexTypeForFields[] accordingly + - Write the c++ class for the new vertex type + + When adding support for new diagram types, do the following: + - Add the new types to diagramTypes + - Write new overloads for CreateDiagramEvaluatorClass[], ContributingDiagramsOfType[] and VerticesForDiagram[] + - Write the necessary c++ code: loop functions, DiagramEvaluator<> specialisations + **********) + +(************* Begin public interface *******************) + +EDMInitialize[] := (subIndexPattern = (Alternatives @@ SARAH`subIndizes[[All, 1]] -> ___);) + +edmFields = Null; +EDMSetEDMFields[fields_List] := (edmFields = fields;) + +EDMCreateFields[] := +Module[{fields, code}, + fields = TreeMasses`GetParticles[]; + + code = (StringJoin @ Riffle[("struct " <> CXXNameOfField[#] <> + ": public Field {\n" <> + TextFormatting`IndentText["static const unsigned numberOfGenerations = " <> + ToString @ TreeMasses`GetDimension[#] <> ";\n"] <> + "};\n" &) /@ fields, "\n"] <> "\n\n" <> + "// Special field families\n" <> + "using Photon = " <> CXXNameOfField @ GetPhoton[] <> ";\n\n" <> + + "// Fields that are their own anti fields\n" <> + StringJoin @ Riffle[("template<> struct " <> + "anti<" <> CXXNameOfField[#] <> ">" <> + " { using type = " <> CXXNameOfField[#] <> "; };" + &) /@ Select[fields, (# == SARAH`AntiField[#] &)], + "\n"] <> "\n\n" <> + + StringJoin @ Riffle[("template<> struct field_indices<" <> + CXXNameOfField[#] <> ">\n" <> + "{\n" <> + IndentText @ + ("using type = std::array + ToString @ Length @ CleanFieldInfo[#][[5]] <> + ">;\n" + ) <> + "};\n" &) /@ fields, "\n"] + ); + + code + ]; + +EDMCreateDiagrams[] := +Module[{code}, + code = StringJoin @ Riffle[(Module[{diagramType = #}, + "template class " <> SymbolName[diagramType] <> ";\n" <> + StringJoin @ Riffle[("template<> class " <> SymbolName[diagramType] <> + "<" <> ToString @ # <> "> {};" + &) /@ diagramSubTypes[diagramType], "\n"] + ] &) /@ diagramTypes, "\n\n"]; + + code = (code <> "\n\n" <> + StringJoin @ Riffle[(Module[{diagramType = #}, + StringJoin @ Riffle[ + ("template\n" <> + "struct DiagramEvaluator<" <> SymbolName[diagramType] <> + "<" <> ToString @ # <> + ">, EDMField, PhotonEmitter, ExchangeField>\n" <> + "{ static double value(const typename field_indices::type &indices, EvaluationContext& context); };" + &) /@ diagramSubTypes[diagramType], "\n\n"]] &) /@ diagramTypes, "\n\n"] + ); + + code + ]; + +EDMCreateInterfaceFunctions[] := +Module[{prototypes, definitions, evaluators}, + evaluators = ConcreteDiagramEvaluators[]; + + prototypes = ("namespace " <> FlexibleSUSY`FSModelName <> "_edm {\n" <> + StringJoin @ Riffle[("double calculate_edm_" <> CXXNameOfField[#] <> + "(" <> + If[TreeMasses`GetDimension[#] =!= 1, + " unsigned generationIndex, ", + " "] <> + "const " <> FlexibleSUSY`FSModelName <> "_mass_eigenstates& model );" + &) /@ edmFields, "\n"] <> + "\n}"); + + definitions = StringJoin @ Riffle[ + Module[{field = #[[1]], fieldEvaluators = #[[2]], + numberOfIndices}, + numberOfIndices = Length @ CleanFieldInfo[field][[5]]; + + "double " <> FlexibleSUSY`FSModelName <> "_edm::calculate_edm_" <> CXXNameOfField[field] <> + "(" <> + If[TreeMasses`GetDimension[field] =!= 1, + " unsigned generationIndex, ", + " "] <> + "const " <> FlexibleSUSY`FSModelName <> "_mass_eigenstates& model )\n" <> + "{\n" <> + IndentText @ + (FlexibleSUSY`FSModelName <> "_mass_eigenstates model_ = model;\n" <> + "EvaluationContext context{ model_ };\n" <> + "std::array + ToString @ numberOfIndices <> + "> indices = {" <> + If[TreeMasses`GetDimension[field] =!= 1, + " generationIndex" <> + 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);" + &) /@ fieldEvaluators, "\n"] <> + "\n\n" <> + "return val;" + ) <> + "\n}"] & /@ evaluators, "\n\n"]; + + {prototypes, definitions} + ]; + +EDMCreateVertexFunctionData[vertexRules_List] := CreateVertices[vertexRules][[1]]; + +EDMCreateDefinitions[vertexRules_List] := +(CreateVertices[vertexRules][[2]] <> "\n\n" <> + CreateEvaluationContextSpecializations[]); + +EDMNPointFunctions[] := +Module[{contributingDiagrams, vertices}, + contributingDiagrams = ContributingDiagrams[]; + + vertices = DeleteDuplicates @ Flatten[VerticesForDiagram /@ + Flatten @ contributingDiagrams[[All, 2]], 1]; + + Flatten[(Null[Null, #] &) /@ ((CouplingsForFields[#] &) /@ vertices)] + ]; + +(**************** End public interface *****************) +Begin["`Private`"]; + +(* The supported vertex types. + They have the same names as their c++ counterparts. *) +vertexTypes = { + SingleComponentedVertex, + LeftAndRightComponentedVertex +}; + +(* The supported diagram types. + They have the same names as their c++ counterparts. *) +diagramTypes = { + OneLoopDiagram +}; + +(* The supported diagram types. + Indexed by the diagram type, gives a set of (c++-compatible) unsigned integer indices. *) +diagramSubTypes[OneLoopDiagram] = { 0, 1 }; (* 0: fermion emits photon, exchange field is a scalar + 1: scalar emits photon, exchange field is a fermion *) + +(**************** CXX conversion routines ***************) + +(* Return a string corresponding to the c++ class name of the field. + Note that "bar" and "conj" get turned into anti<...>::type! *) +CXXNameOfField[p_] := SymbolName[p]; +CXXNameOfField[SARAH`bar[p_]] := "anti<" <> SymbolName[p] <> ">::type"; +CXXNameOfField[Susyno`LieGroups`conj[p_]] := "anti<" <> SymbolName[p] <> ">::type"; + +(**************** Other Functions ***************) + +GetPhoton[] := SARAH`Photon; + +IsLorentzIndex[index_] := StringMatchQ[ToString @ index, "lt" ~~ __]; + +StripLorentzIndices[p_Symbol] := p; +StripLorentzIndices[SARAH`bar[p_]] := SARAH`bar[StripLorentzIndices[p]]; +StripLorentzIndices[Susyno`LieGroups`conj[p_]] := Susyno`LieGroups`conj[StripLorentzIndices[p]]; +StripLorentzIndices[p_] := Module[{remainingIndices}, + remainingIndices = Select[p[[1]], (!IsLorentzIndex[#] &)]; + If[Length[remainingIndices] === 0, Head[p], + Head[p][remainingIndices]] + ]; + +CleanFieldInfo[field_] := Module[{fieldInfo = Cases[SARAH`Particles[FlexibleSUSY`FSEigenstates], + {SARAH`getParticleName @ field, ___}][[1]]}, + fieldInfo = DeleteCases[fieldInfo, {SARAH`generation, 1}, {2}]; + DeleteCases[fieldInfo, {SARAH`lorentz, _}, {2}] + ]; + +CreateEvaluationContextSpecializations[] := +Module[{fields, contributingDiagrams, photonEmitters}, + fields = Select[TreeMasses`GetParticles[], (! TreeMasses`IsGhost[#] &)]; + + contributingDiagrams = ContributingDiagrams[]; + photonEmitters = DeleteDuplicates @ Flatten[contributingDiagrams[[All, 2]], 1][[All, 3]]; + + StringJoin @ Riffle[(Module[{fieldInfo = CleanFieldInfo[#], numberOfIndices}, + numberOfIndices = Length @ fieldInfo[[5]]; + + "template<> double EvaluationContext::mass<" <> ToString[#] <> + ">( const std::array ToString @ numberOfIndices <> + "> &indices ) const\n" <> + "{ return model.get_M" <> CXXNameOfField[#] <> + If[TreeMasses`GetDimension[#] === 1, "()", "( indices[0] )"] <> "; }" + ] &) /@ fields, "\n\n"] <> "\n\n" <> + + StringJoin @ Riffle[(Module[{fieldInfo = CleanFieldInfo[#], + photonVertexType = VertexTypeForFields[{GetPhoton[], #, SARAH`AntiField @ #}], + numberOfIndices}, + numberOfIndices = Length @ fieldInfo[[5]]; + + "template<>\n" <> + "double EvaluationContext::charge<" <> CXXNameOfField[#] <> + ">( const std::array ToString @ numberOfIndices <> + "> &indices ) const\n" <> + "{\n" <> + IndentText @ + ("using PhotonVertex = VertexFunction + CXXNameOfField[#] <> ", " <> CXXNameOfField[SARAH`AntiField @ #] <> + ">;\n\n" <> + "return PhotonVertex::vertex( concatenate( indices, indices ), *this )" <> + If[photonVertexType === SingleComponentedVertex, + ".value().real();\n", + ".left().real();\n"] + ) <> + "}"] &) /@ photonEmitters, "\n\n"] + ]; + +(* Find all diagrams of the type type_, testing all corresponding combinations of fields *) +(* IMPORTANT: Return value should have the format + {{edmField1, {Diagram[DIAGRAMTYPENAME[_Integer], Fields___], Diagram[...], ...}}, + {edmField2, {...}}, + ...} *) +ContributingDiagramsOfType[OneLoopDiagram] := + Module[{edmField = #, diagrams = SARAH`InsFields[ + {{C[#, SARAH`AntiField[SARAH`FieldToInsert[1]], + SARAH`AntiField[SARAH`FieldToInsert[2]]], + C[SARAH`FieldToInsert[1], GetPhoton[], + SARAH`AntiField[SARAH`FieldToInsert[1]]], + C[SARAH`FieldToInsert[1], SARAH`FieldToInsert[2], + SARAH`AntiField[#]]}, + {SARAH`FieldToInsert[1], SARAH`FieldToInsert[2]}}], + subtypedDiagrams, uniqueDiagrams}, + + If[TreeMasses`IsFermion[edmField] =!= True, + {edmField,{}}, + + subtypedDiagrams = (Module[{photonEmitter = #[[2,1]], + exchangeField = #[[2,2]], + subType}, + subType = If[TreeMasses`IsFermion[photonEmitter] && + TreeMasses`IsScalar[exchangeField], + 0, + If[TreeMasses`IsScalar[photonEmitter] && + TreeMasses`IsFermion[exchangeField], + 1]]; + If[subType === Null, + Null, + Diagram[OneLoopDiagram[subType], edmField, photonEmitter, exchangeField]] + ] + &) /@ diagrams; + uniqueDiagrams = DeleteDuplicates @ Cases[subtypedDiagrams, Except[Null]]; + {edmField, uniqueDiagrams}]] & /@ edmFields; + +(* Returns the necessary c++ code corresponding to the vertices that need to be calculated. + The returned value is a list {prototypes, definitions}. *) +CreateVertices[vertexRules_List] := + Module[{contributingDiagrams, vertices, + vertexClassesPrototypes, vertexClassesDefinitions}, + contributingDiagrams = ContributingDiagrams[]; + + vertices = DeleteDuplicates @ Flatten[VerticesForDiagram /@ + Flatten @ contributingDiagrams[[All, 2]], 1]; + + If[vertices === {}, + {"",""}, + + {vertexClassesPrototypes, vertexClassesDefinitions} = Transpose @ + ((CreateVertexFunction[#, vertexRules] &) /@ vertices); + (StringJoin @ Riffle[#, "\n\n"] &) /@ {vertexClassesPrototypes, vertexClassesDefinitions} + ]]; + +(* Returns the vertices that are present in the specified diagram. + This function should be overloaded for future diagram types. *) +VerticesForDiagram[Diagram[loopDiagram_OneLoopDiagram, edmField_, photonEmitter_, exchangeField_]] := + Module[{edmVertex, photonVertex}, + edmVertex = {SARAH`AntiField[edmField], photonEmitter, exchangeField}; + photonVertex = {GetPhoton[], photonEmitter, SARAH`AntiField[photonEmitter]}; + {edmVertex, photonVertex} + ]; + +(* Returns the vertex type for a vertex with a given list of fields *) +VertexTypeForFields[fields_List] := + Module[{fermions, scalarCount, vectorCount, fermionCount, vertexType = Null}, + fermions = Select[fields, TreeMasses`IsFermion]; + + scalarCount = Length @ Select[fields, TreeMasses`IsScalar]; + vectorCount = Length @ Select[fields, TreeMasses`IsVector]; + fermionCount = Length @ fermions; + + If[fermionCount === 2 && scalarCount === 1 && vectorCount === 0, + vertexType = LeftAndRightComponentedVertex]; + If[fermionCount === 2 && scalarCount === 0 && vectorCount === 1, + If[fermions[[1]] === SARAH`AntiField[fermions[[2]]], + vertexType = LeftAndRightComponentedVertex]]; + If[fermionCount === 0 && scalarCount === 2 && vectorCount === 1, + vertexType = SingleComponentedVertex]; + + vertexType + ]; + +(* Returns the different SARAH`Cp coupling parts for a vertex with a given list of fields *) +CouplingsForFields[fields_List] := + Module[{vertexType, couplings}, + vertexType = VertexTypeForFields[fields]; + couplings = {SARAH`Cp @@ fields}; + + Switch[vertexType, + SingleComponentedVertex, couplings, + LeftAndRightComponentedVertex, {couplings[[1]][SARAH`PL], couplings[[1]][SARAH`PR]}] + ]; + +(* Creates the actual c++ code for a vertex with given fields. + This involves creating the VertexFunctionData<> code as well as + the VertexFunction<> code. You should never need to change this code! *) +CreateVertexFunction[fields_List, vertexRules_List] := + Module[{prototype, definition, + parsedVertex, dataClassName, functionClassName, fieldIndexStartF, + fieldIndexStart, indexBounds}, + parsedVertex = ParseVertex[fields, vertexRules]; + + dataClassName = "VertexFunctionData<" <> StringJoin @ Riffle[CXXNameOfField /@ fields, ", "] <> ">"; + functionClassName = "VertexFunction<" <> StringJoin @ Riffle[CXXNameOfField /@ fields, ", "] <> ">"; + + fieldIndexStartF[1] = 0; + fieldIndexStartF[pIndex_] := fieldIndexStartF[pIndex-1] + NumberOfIndices[parsedVertex, pIndex-1]; + fieldIndexStartF[Length[fields]+1] = NumberOfIndices[parsedVertex]; + + fieldIndexStart = Table[fieldIndexStartF[i], {i, 1, Length[fields] + 1}]; + + indexBounds = IndexBounds[parsedVertex]; + + prototype = ("template<> struct " <> dataClassName <> "\n" <> + "{\n" <> + IndentText @ + ("static constexpr IndexBounds<" <> ToString @ NumberOfIndices[parsedVertex] <> + "> index_bounds" <> + If[NumberOfIndices[parsedVertex] =!= 0, + " = { " <> + "{ " <> StringJoin @ Riffle[ToString /@ indexBounds[[1]], ", "] <> " }, " <> + "{ " <> StringJoin @ Riffle[ToString /@ indexBounds[[2]], ", "] <> " } };\n" + , + "{};\n" + ] <> + "static constexpr unsigned fieldIndexStart[" <> ToString @ Length[fieldIndexStart] <> + "] = { " <> StringJoin @ Riffle[ToString /@ fieldIndexStart, ", "] <> + " };\n" <> + "using vertex_type = " <> VertexClassName[parsedVertex] <> ";\n" + ) <> + "};"); + + definition = ("template<> " <> functionClassName <> "::vertex_type\n" <> + functionClassName <> "::vertex(const indices_type &indices, const EvaluationContext &context)\n" <> + "{\n" <> + IndentText @ VertexFunctionBody[parsedVertex] <> "\n" <> + "}"); + + {prototype, definition} + ]; + +(* Creates local declarations of field indices, whose values are taken + from the elements of `arrayName'. + *) +DeclareIndices[indexedFields_List, arrayName_String] := + Module[{p, total = 0, fieldIndexList, decl = "", idx}, + DeclareIndex[idx_, num_Integer, an_String] := ( + "const unsigned " <> CConversion`ToValidCSymbolString[idx] <> + " = " <> an <> "[" <> ToString[num] <> "];\n"); + For[p = 1, p <= Length[indexedFields], p++, + fieldIndexList = Vertices`FieldIndexList[indexedFields[[p]]]; + decl = decl <> StringJoin[DeclareIndex[#, total++, arrayName]& /@ fieldIndexList]; + ]; + Assert[total == Total[Length[Vertices`FieldIndexList[#]]& /@ indexedFields]]; + decl + ]; + +GetComplexScalarCType[] := + CConversion`CreateCType[CConversion`ScalarType[CConversion`complexScalarCType]]; + +(* ParsedVertex structure: + ParsedVertex[ + {numP1Indices, numP2Indices, ...}, + {{minIndex1, minIndex2, ...}, {maxIndex1+1, maxIndex2+1, ...}}, + VertexClassName, + VertexFunctionBody + ] + + Getters are available! Given below ParseVertex[] + *) + +(* The heart of the algorithm! From the field content, determine all + necessary information. *) +ParseVertex[fields_List, vertexRules_List] := + Module[{indexedFields, numberOfIndices, declareIndices, + parsedVertex, vertexClassName, vertexFunctionBody, + fieldInfo, trIndexBounds, indexBounds, + expr, exprL, exprR}, + indexedFields = MapIndexed[(Module[{field = #1, + index = #2[[1]]}, + SARAH`getFull[#1] /. SARAH`subGC[index] /. SARAH`subIndFinal[index,index] + ] &), fields]; + indexedFields = StripLorentzIndices /@ indexedFields; + + + numberOfIndices = ((Length @ Vertices`FieldIndexList[#] &) /@ indexedFields); + declareIndices = DeclareIndices[indexedFields, "indices"]; + + vertexClassName = SymbolName[VertexTypeForFields[fields]]; + vertexFunctionBody = Switch[vertexClassName, + "SingleComponentedVertex", + expr = (SARAH`Cp @@ fields) /. vertexRules; + expr = TreeMasses`ReplaceDependenciesReverse[expr]; + declareIndices <> + Parameters`CreateLocalConstRefs[expr] <> "\n" <> + "const " <> GetComplexScalarCType[] <> " result = " <> + Parameters`ExpressionToString[expr] <> ";\n\n" <> + "return vertex_type(result);", + + "LeftAndRightComponentedVertex", + exprL = SARAH`Cp[Sequence @@ fields][SARAH`PL] /. vertexRules; + exprR = SARAH`Cp[Sequence @@ fields][SARAH`PR] /. vertexRules; + exprL = TreeMasses`ReplaceDependenciesReverse[exprL]; + exprR = TreeMasses`ReplaceDependenciesReverse[exprR]; + declareIndices <> + Parameters`CreateLocalConstRefs[exprL + exprR] <> "\n" <> + "const " <> GetComplexScalarCType[] <> " left = " <> + Parameters`ExpressionToString[exprL] <> ";\n\n" <> + "const " <> GetComplexScalarCType[] <> " right = " <> + Parameters`ExpressionToString[exprR] <> ";\n\n" <> + "return vertex_type(left, right);"]; + + fieldInfo = CleanFieldInfo /@ fields; + + trIndexBounds = Cases[Flatten[(With[{fieldIndex = #}, + (If[#[[1]] === SARAH`generation, + {fieldInfo[[fieldIndex, 2]]-1, fieldInfo[[fieldIndex, 3]]}, + {1, #[[2]]}] + &) /@ fieldInfo[[fieldIndex, 5]]] + &) /@ Table[i, {i, Length[fields]}], + 1], + Except[{}]]; + + If[trIndexBounds === {}, + indexBounds = {{},{}}, + indexBounds = Transpose @ trIndexBounds]; + + parsedVertex = ParsedVertex[numberOfIndices, + indexBounds, + vertexClassName, + vertexFunctionBody]; + + parsedVertex + ]; + +(** Getters to the ParsedVertex structure **) +NumberOfIndices[parsedVertex_ParsedVertex] := Total @ parsedVertex[[1]]; +NumberOfIndices[parsedVertex_ParsedVertex, pIndex_Integer] := parsedVertex[[1, pIndex]]; + +IndexBounds[parsedVertex_ParsedVertex] := parsedVertex[[2]]; + +VertexClassName[parsedVertex_ParsedVertex] := parsedVertex[[3]]; +VertexFunctionBody[parsedVertex_ParsedVertex] := parsedVertex[[4]]; +(** End getters **) + +(* Find all contributing diagrams *) +cachedContributingDiagrams = Null; +ContributingDiagrams[] := + Module[{diagrams}, + If[cachedContributingDiagrams =!= Null, + cachedContributingDiagrams, + + LoadVerticesIfNecessary[]; + + diagrams = Flatten[(ContributingDiagramsOfType[#] &) + /@ diagramTypes + , 1]; + cachedContributingDiagrams = ({#, Union @ + (Sequence @@ Cases[diagrams, + {#, diags_List} -> diags])} &) /@ edmFields; + + cachedContributingDiagrams] + ]; + +LoadVerticesIfNecessary[] := + Module[{}, + If[SARAH`VertexList3 =!= List || Length[SARAH`VertexList3] === 0, + SA`CurrentStates = FlexibleSUSY`FSEigenstates; + SARAH`InitVertexCalculation[FlexibleSUSY`FSEigenstates, False]; + SARAH`ReadVertexList[FlexibleSUSY`FSEigenstates, False, False, True]; + SARAH`MakeCouplingLists; + ]; + ]; + +ConcreteDiagramEvaluators[] := + ({#[[1]], + (("DiagramEvaluator<" <> SymbolName @ Head @ #[[1]] <> "<" <> + ToString @ #[[1,1]] <> ">, " <> + StringJoin @ (Riffle[CXXNameOfField /@ List @@ #[[2;;]], ", "]) <> + ">" &) + /@ #[[2]]) } &) /@ ContributingDiagrams[]; + +End[]; + +EndPackage[]; diff --git a/meta/FSMathLink.m b/meta/FSMathLink.m index 3235d3a2a..85ca0b504 100644 --- a/meta/FSMathLink.m +++ b/meta/FSMathLink.m @@ -164,9 +164,20 @@ PutSpectrum[pars_List, link_String] := StringJoin[PutParameter[#,link]& /@ pars]; -PutObservable[obs_, type_, link_String] := - "MLPutRuleTo(" <> link <> ", OBSERVABLE(" <> Observables`GetObservableName[obs] <> - "), \"" <> ToString[obs] <> "\");\n"; +ObsToStr[obs_?NumericQ] := ToString[obs]; +ObsToStr[obs_] := "\"" <> ToString[obs] <> "\""; + +HeadToStr[sym_] := "\"" <> ToString[sym] <> "\""; +HeadsToStr[{}] := ""; +HeadsToStr[l_List] := ", {" <> StringJoin[Riffle[HeadToStr /@ l, ", "]] <> "}"; + +PutObservable[obs_[sub_], type_, link_String, heads_:{}] := + PutObservable[sub, type, link, Join[heads, {obs}]]; + +PutObservable[obs_, type_, link_String, heads_:{}] := + "MLPutRuleTo(" <> link <> ", OBSERVABLE(" <> + Observables`GetObservableName[Composition[Sequence @@ heads][obs]] <> + "), " <> ObsToStr[obs] <> HeadsToStr[heads] <> ");\n"; PutObservables[obs_List, link_String] := StringJoin[PutObservable[#, Observables`GetObservableType[#], link]& /@ obs]; diff --git a/meta/FlexibleSUSY.m b/meta/FlexibleSUSY.m index 80c888fc7..2fa838400 100644 --- a/meta/FlexibleSUSY.m +++ b/meta/FlexibleSUSY.m @@ -23,6 +23,7 @@ "ThreeLoopMSSM`", "Observables`", "GMuonMinus2`", + "EDM`", "EffectiveCouplings`", "FlexibleEFTHiggsMatching`", "FSMathLink`", @@ -1756,12 +1757,12 @@ corresponding tadpole is real or imaginary (only in models with CP WriteGMuonMinus2Class[vertexRules_List, files_List] := Module[{particles, muonFunctionPrototypes, diagrams, vertexFunctionData, definitions, calculationCode, getMSUSY, getQED2L}, - particles = GMuonMinus2`CreateParticles[]; - muonFunctionPrototypes = GMuonMinus2`CreateMuonFunctions[vertexRules][[1]]; - diagrams = GMuonMinus2`CreateDiagrams[]; - vertexFunctionData = GMuonMinus2`CreateVertexFunctionData[vertexRules]; - definitions = GMuonMinus2`CreateDefinitions[vertexRules]; - calculationCode = GMuonMinus2`CreateCalculation[]; + particles = GMuonMinus2`GMuonMinus2CreateParticles[]; + muonFunctionPrototypes = GMuonMinus2`GMuonMinus2CreateMuonFunctions[vertexRules][[1]]; + diagrams = GMuonMinus2`GMuonMinus2CreateDiagrams[]; + vertexFunctionData = GMuonMinus2`GMuonMinus2CreateVertexFunctionData[vertexRules]; + definitions = GMuonMinus2`GMuonMinus2CreateDefinitions[vertexRules]; + calculationCode = GMuonMinus2`GMuonMinus2CreateCalculation[]; getMSUSY = GMuonMinus2`GetMSUSY[]; getQED2L = GMuonMinus2`GetQED2L[]; @@ -1778,6 +1779,26 @@ corresponding tadpole is real or imaginary (only in models with CP } ]; ]; +(* Write the EDM c++ files *) +WriteEDMClass[vertexRules_List, files_List] := + Module[{fields, diagrams, interfaceFunctions, + vertexFunctionData, definitions}, + fields = EDM`EDMCreateFields[]; + diagrams = EDM`EDMCreateDiagrams[]; + interfaceFunctions = EDM`EDMCreateInterfaceFunctions[]; + vertexFunctionData = EDM`EDMCreateVertexFunctionData[vertexRules]; + definitions = EDM`EDMCreateDefinitions[vertexRules]; + WriteOut`ReplaceInFiles[files, + { "@EDM_Fields@" -> fields, + "@EDM_Diagrams@" -> diagrams, + "@EDM_VertexFunctionData@" -> vertexFunctionData, + "@EDM_Definitions@" -> definitions, + "@EDM_InterfaceFunctionPrototypes@" -> interfaceFunctions[[1]], + "@EDM_InterfaceFunctionDefinitions@" -> interfaceFunctions[[2]], + Sequence @@ GeneralReplacementRules[] + } ]; + ]; + GetBVPSolverHeaderName[solver_] := Switch[solver, FlexibleSUSY`TwoScaleSolver, "two_scale", @@ -2292,7 +2313,7 @@ corresponding tadpole is real or imaginary (only in models with CP ]; (* Get all nPointFunctions that GMM2 needs *) -PrepareGMuonMinus2[] := GMuonMinus2`NPointFunctions[]; +PrepareGMuonMinus2[] := GMuonMinus2`GMuonMinus2NPointFunctions[]; PrepareUnrotatedParticles[eigenstates_] := Module[{nonMixedParticles = {}, nonMixedParticlesFile}, @@ -2724,7 +2745,7 @@ corresponding tadpole is real or imaginary (only in models with CP MakeFlexibleSUSY[OptionsPattern[]] := Module[{nPointFunctions, runInputFile, initialGuesserInputFile, - gmm2Vertices = {}, + gmm2Vertices = {}, edmFields, susyBetaFunctions, susyBreakingBetaFunctions, numberOfSusyParameters, anomDim, inputParameters (* list of 3-component lists of the form {name, block, type} *), @@ -3579,7 +3600,19 @@ corresponding tadpole is real or imaginary (only in models with CP {{FileNameJoin[{$flexiblesusyTemplateDir, "a_muon.hpp.in"}], FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> "_a_muon.hpp"}]}, {FileNameJoin[{$flexiblesusyTemplateDir, "a_muon.cpp.in"}], - FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> "_a_muon.cpp"}]}}]; + FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> "_a_muon.cpp"}]}}]; + + Print["Creating class EDM"]; + edmFields = DeleteDuplicates @ Cases[Observables`GetRequestedObservables[extraSLHAOutputBlocks], + FlexibleSUSYObservable`EDM[p_[__]|p_] :> p]; + EDM`EDMInitialize[]; + EDM`EDMSetEDMFields[edmFields]; + + WriteEDMClass[Join[vertexRules, Vertices`VertexRules[EDM`EDMNPointFunctions[], Lat$massMatrices]], + {{FileNameJoin[{$flexiblesusyTemplateDir, "edm.hpp.in"}], + FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> "_edm.hpp"}]}, + {FileNameJoin[{$flexiblesusyTemplateDir, "edm.cpp.in"}], + FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> "_edm.cpp"}]}}]; PrintHeadline["Creating Mathematica interface"]; Print["Creating LibraryLink ", FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> ".mx"}], " ..."]; diff --git a/meta/GMuonMinus2.m b/meta/GMuonMinus2.m index 278ec11a7..eceacc448 100644 --- a/meta/GMuonMinus2.m +++ b/meta/GMuonMinus2.m @@ -1,15 +1,15 @@ BeginPackage["GMuonMinus2`", {"SARAH`", "CConversion`", "TextFormatting`", "TreeMasses`", "LoopMasses`", "Vertices`"}]; -CreateParticles::usage="Returns the c++ code that contains all particle classes"; -CreateMuonFunctions::usage="Returns the c++ code that contains all muon functions"; -CreateDiagrams::usage="Returns the c++ code that contains all relevant diagram classes"; -CreateVertexFunctionData::usage="Returns the c++ code that contains all relevant vertex function data"; +GMuonMinus2CreateParticles::usage="Returns the c++ code that contains all particle classes"; +GMuonMinus2CreateMuonFunctions::usage="Returns the c++ code that contains all muon functions"; +GMuonMinus2CreateDiagrams::usage="Returns the c++ code that contains all relevant diagram classes"; +GMuonMinus2CreateVertexFunctionData::usage="Returns the c++ code that contains all relevant vertex function data"; -CreateCalculation::usage="Returns the c++ code that performs the actual calculation the magnetic moment"; +GMuonMinus2CreateCalculation::usage="Returns the c++ code that performs the actual calculation the magnetic moment"; -CreateDefinitions::usage="Returns the c++ that contains all function definitions" +GMuonMinus2CreateDefinitions::usage="Returns the c++ that contains all function definitions" -NPointFunctions::usage="Returns a list of all n point functions that are needed. Actually it is a list of fake functions to extract vertex functions..."; +GMuonMinus2NPointFunctions::usage="Returns a list of all n point functions that are needed. Actually it is a list of fake functions to extract vertex functions..."; GetMSUSY::usage="returns minimum charged BSM mass"; @@ -49,7 +49,7 @@ If you add new kinds of vertices (e.g for new diagram types): ]; (* Create c++ classes for all particles *) -CreateParticles[] := +GMuonMinus2CreateParticles[] := Module[{particles, code}, (* Get a list of all particles *) particles = TreeMasses`GetParticles[]; @@ -111,7 +111,7 @@ If you add new kinds of vertices (e.g for new diagram types): ]; muonFunctions = Null; -CreateMuonFunctions[vertexRules_List] := +GMuonMinus2CreateMuonFunctions[vertexRules_List] := Module[{muonIndex, muonFamily, prototypes = "", definitions, contextMuonPole}, If[muonFunctions =!= Null, Return[muonFunctions]]; @@ -145,7 +145,7 @@ If you add new kinds of vertices (e.g for new diagram types): {prototypes, definitions} ]; -CreateDiagrams[] := +GMuonMinus2CreateDiagrams[] := Module[{diagramTypes, diagramTypeHeads, code}, diagrams = contributingFeynmanDiagramTypes; diagramHeads = DeleteDuplicates @ (Head /@ diagrams); @@ -167,7 +167,7 @@ If you add new kinds of vertices (e.g for new diagram types): Return[code]; ]; -CreateVertexFunctionData[vertexRules_List] := CreateVertices[vertexRules][[1]]; +GMuonMinus2CreateVertexFunctionData[vertexRules_List] := CreateVertices[vertexRules][[1]]; CreateDiagramEvaluatorClass[type_OneLoopDiagram] := ("template\n" <> @@ -206,7 +206,7 @@ If you add new kinds of vertices (e.g for new diagram types): "return qed_2L;"; calculationCode = Null; -CreateCalculation[] := +GMuonMinus2CreateCalculation[] := Module[{code}, (* If we have been here before return the old result *) If[calculationCode =!= Null, Return[calculationCode]]; @@ -228,13 +228,13 @@ If you add new kinds of vertices (e.g for new diagram types): Return[code]; ]; -CreateDefinitions[vertexRules_List] := +GMuonMinus2CreateDefinitions[vertexRules_List] := (CreateEvaluationContextSpecializations[] <> "\n\n" <> - CreateMuonFunctions[vertexRules][[2]] <> "\n\n" <> + GMuonMinus2CreateMuonFunctions[vertexRules][[2]] <> "\n\n" <> CreateVertices[vertexRules][[2]]); nPointFunctions = Null; -NPointFunctions[] := +GMuonMinus2NPointFunctions[] := Module[{contributingDiagrams, vertices}, If[nPointFunctions =!= Null, Return[nPointFunctions]]; diff --git a/meta/Observables.m b/meta/Observables.m index cfcabb879..c1e3c414e 100644 --- a/meta/Observables.m +++ b/meta/Observables.m @@ -74,8 +74,8 @@ GetObservableDescription[obs_ /; obs === FlexibleSUSYObservable`CpHiggsGluonGluon] := "effective H-Gluon-Gluon coupling"; GetObservableDescription[obs_ /; obs === FlexibleSUSYObservable`CpPseudoScalarPhotonPhoton] := "effective A-Photon-Photon coupling"; GetObservableDescription[obs_ /; obs === FlexibleSUSYObservable`CpPseudoScalarGluonGluon] := "effective A-Gluon-Gluon coupling"; -GetObservableDescription[FlexibleSUSYObservable`EDM[p_[idx_]]] := GetObservableDescription[FlexibleSUSYObservable`EDM[p]] <> "(" <> ToString[idx] <> ")"; -GetObservableDescription[FlexibleSUSYObservable`EDM[p_]] := "electric dipole moment of " <> CConversion`ToValidCSymbolString[p]; +GetObservableDescription[FlexibleSUSYObservable`EDM[p_[idx_]]] := "electric dipole moment of " <> CConversion`ToValidCSymbolString[p] <> "(" <> ToString[idx] <> ") [1/GeV]"; +GetObservableDescription[FlexibleSUSYObservable`EDM[p_]] := "electric dipole moment of " <> CConversion`ToValidCSymbolString[p] <> " [1/GeV]"; GetObservableType[obs_ /; obs === FlexibleSUSYObservable`aMuon] := CConversion`ScalarType[CConversion`realScalarCType]; GetObservableType[obs_ /; obs === FlexibleSUSYObservable`aMuonUncertainty] := CConversion`ScalarType[CConversion`realScalarCType]; diff --git a/meta/WriteOut.m b/meta/WriteOut.m index 2e1cfeb5d..414cc0ced 100644 --- a/meta/WriteOut.m +++ b/meta/WriteOut.m @@ -885,13 +885,13 @@ CConversion`ToValidCSymbolString[c] <> "_slha"; GetSLHATrilinearCouplingType[c_] := - CConversion`ToRealType[Parameters`GetType[c]]; + Parameters`GetType[c]; CreateSLHASoftSquaredMassName[c_] := CConversion`ToValidCSymbolString[c] <> "_slha"; GetSLHASoftSquaredMassType[c_] := - CConversion`ToRealType[Parameters`GetType[c]]; + Parameters`GetType[c]; CreateSLHAYukawaDefinition[] := StringJoin[Parameters`CreateParameterDefinitionAndDefaultInitialize[ @@ -976,19 +976,24 @@ If[Parameters`IsOutputParameter[{vL, vR}] && ParametersHaveSameDimension[{vL, vR, #}], result = result <> - CreateSLHATrilinearCouplingName[#] <> " = (" <> - CreateSLHAFermionMixingMatrixName[vR] <> ".conjugate() * " <> - "MODELPARAMETER(" <> - CConversion`ToValidCSymbolString[#] <> ") * " <> - CreateSLHAFermionMixingMatrixName[vL] <> ".adjoint()" <> - ").real();\n"; + CreateSLHATrilinearCouplingName[#] <> " = " <> + CConversion`CastTo[ + CreateSLHAFermionMixingMatrixName[vR] <> ".conjugate() * " <> + "MODELPARAMETER(" <> + CConversion`ToValidCSymbolString[#] <> ") * " <> + CreateSLHAFermionMixingMatrixName[vL] <> ".adjoint()", + GetSLHATrilinearCouplingType[#] + ] <> ";\n"; , Print["Warning: Cannot convert Trilinear coupling ", #, " to SLHA, because ", {vL,vR}, " are not defined", " or have incompatible dimension."]; result = result <> - CreateSLHATrilinearCouplingName[#] <> " = MODELPARAMETER(" <> - CConversion`ToValidCSymbolString[#] <> ").real();\n"; + CreateSLHATrilinearCouplingName[#] <> " = " <> + CConversion`CastTo[ + "MODELPARAMETER(" <> CConversion`ToValidCSymbolString[#] <> ")", + GetSLHATrilinearCouplingType[#] + ] <> ";\n"; ]; ]& /@ tril; result @@ -1038,28 +1043,35 @@ ParametersHaveSameDimension[{vL, vR, #}], If[IsLeftHanded[#], result = result <> - CreateSLHASoftSquaredMassName[#] <> " = (" <> - CreateSLHAFermionMixingMatrixName[vL] <> " * " <> - "MODELPARAMETER(" <> - CConversion`ToValidCSymbolString[#] <> ") * " <> - CreateSLHAFermionMixingMatrixName[vL] <> ".adjoint()" <> - ").real();\n"; + CreateSLHASoftSquaredMassName[#] <> " = " <> + CConversion`CastTo[ + CreateSLHAFermionMixingMatrixName[vL] <> " * " <> + "MODELPARAMETER(" <> + CConversion`ToValidCSymbolString[#] <> ") * " <> + CreateSLHAFermionMixingMatrixName[vL] <> ".adjoint()", + GetSLHASoftSquaredMassType[#] + ] <> ";\n"; , result = result <> - CreateSLHASoftSquaredMassName[#] <> " = (" <> - CreateSLHAFermionMixingMatrixName[vR] <> ".conjugate() * " <> - "MODELPARAMETER(" <> - CConversion`ToValidCSymbolString[#] <> ") * " <> - CreateSLHAFermionMixingMatrixName[vR] <> ".transpose()" <> - ").real();\n"; + CreateSLHASoftSquaredMassName[#] <> " = " <> + CConversion`CastTo[ + CreateSLHAFermionMixingMatrixName[vR] <> ".conjugate() * " <> + "MODELPARAMETER(" <> + CConversion`ToValidCSymbolString[#] <> ") * " <> + CreateSLHAFermionMixingMatrixName[vR] <> ".transpose()", + GetSLHASoftSquaredMassType[#] + ] <> ";\n"; ]; , Print["Warning: Cannot convert soft squared mass ", #, " to SLHA, because ", {vL,vR}, " are not defined", " or have incompatible dimension."]; result = result <> - CreateSLHASoftSquaredMassName[#] <> " = MODELPARAMETER(" <> - CConversion`ToValidCSymbolString[#] <> ").real();\n"; + CreateSLHASoftSquaredMassName[#] <> " = " <> + CConversion`CastTo[ + "MODELPARAMETER(" <> CConversion`ToValidCSymbolString[#] <> ")", + GetSLHASoftSquaredMassType[#] + ] <> ";\n"; ]; ]& /@ massSq; result diff --git a/model_files/CMSSMCPV/FlexibleSUSY.m.in b/model_files/CMSSMCPV/FlexibleSUSY.m.in index 6bb6566c7..6053466f4 100644 --- a/model_files/CMSSMCPV/FlexibleSUSY.m.in +++ b/model_files/CMSSMCPV/FlexibleSUSY.m.in @@ -85,7 +85,10 @@ ExtraSLHAOutputBlocks = { {1, Hold[SUSYScale]}, {2, Hold[LowScale]} } }, {FlexibleSUSYLowEnergy, - {{21, FlexibleSUSYObservable`aMuon} } }, + {{21, FlexibleSUSYObservable`aMuon}, + {23, FlexibleSUSYObservable`EDM[Fe[1]]}, + {24, FlexibleSUSYObservable`EDM[Fe[2]]}, + {25, FlexibleSUSYObservable`EDM[Fe[3]]} } }, {EFFHIGGSCOUPLINGS, NoScale, {{1, FlexibleSUSYObservable`CpHiggsPhotonPhoton}, {2, FlexibleSUSYObservable`CpHiggsGluonGluon}, diff --git a/model_files/MSSMCPV/FlexibleSUSY.m.in b/model_files/MSSMCPV/FlexibleSUSY.m.in index 734bfca2e..849c4b8be 100644 --- a/model_files/MSSMCPV/FlexibleSUSY.m.in +++ b/model_files/MSSMCPV/FlexibleSUSY.m.in @@ -126,7 +126,10 @@ ExtraSLHAOutputBlocks = { {1, Hold[SUSYScale]}, {2, Hold[LowScale]} } }, {FlexibleSUSYLowEnergy, - {{21, FlexibleSUSYObservable`aMuon} } }, + {{21, FlexibleSUSYObservable`aMuon}, + {23, FlexibleSUSYObservable`EDM[Fe[1]]}, + {24, FlexibleSUSYObservable`EDM[Fe[2]]}, + {25, FlexibleSUSYObservable`EDM[Fe[3]]} } }, {EFFHIGGSCOUPLINGS, NoScale, {{1, FlexibleSUSYObservable`CpHiggsPhotonPhoton}, {2, FlexibleSUSYObservable`CpHiggsGluonGluon}, diff --git a/model_files/NMSSMCPV/FlexibleSUSY.m.in b/model_files/NMSSMCPV/FlexibleSUSY.m.in index 667e2d19b..6adc274a9 100644 --- a/model_files/NMSSMCPV/FlexibleSUSY.m.in +++ b/model_files/NMSSMCPV/FlexibleSUSY.m.in @@ -98,6 +98,11 @@ ExtraSLHAOutputBlocks = { {{0, Hold[HighScale]}, {1, Hold[SUSYScale]}, {2, Hold[LowScale]} } }, + {FlexibleSUSYLowEnergy, + {{21, FlexibleSUSYObservable`aMuon}, + {23, FlexibleSUSYObservable`EDM[Fe[1]]}, + {24, FlexibleSUSYObservable`EDM[Fe[2]]}, + {25, FlexibleSUSYObservable`EDM[Fe[3]]} } }, {NMSSMRUN, {{1 , Re[\[Lambda]] }, {2 , Re[\[Kappa]] }, diff --git a/src/ew_input.hpp b/src/ew_input.hpp index da28676dc..75f1487af 100644 --- a/src/ew_input.hpp +++ b/src/ew_input.hpp @@ -80,7 +80,7 @@ namespace Electroweak_constants { const double PMNS_DELTA = 0.; const double PMNS_ALPHA1 = 0.; const double PMNS_ALPHA2 = 0.; - const double gfermi = 1.16637e-5; ///< Fermi constant G_F + const double gfermi = 1.1663787e-5; ///< Fermi constant G_F const double yeSM = 2.85784368E-06; ///< Ye(1,1) MS-bar in the SM at Q = MZ const double ymSM = 5.90911374E-04; ///< Ye(2,2) MS-bar in the SM at Q = MZ const double ylSM = 9.95869693E-03; ///< Ye(3,3) MS-bar in the SM at Q = MZ diff --git a/src/mathlink_utils.hpp b/src/mathlink_utils.hpp index 0ba023145..356803f86 100644 --- a/src/mathlink_utils.hpp +++ b/src/mathlink_utils.hpp @@ -125,8 +125,22 @@ inline void MLPutRule(MLINK link, const std::string& name, const std::vector(name.c_str()), name.size()); } -template -void MLPutRuleTo(MLINK link, T t, const std::string& name, const std::vector& heads = {}) +inline void MLPutRule(MLINK link, int number, const std::vector& heads = {}) +{ + MLPutFunction(link, "Rule", 2); + MLPutHeads(link, heads); + MLPutInteger(link, number); +} + +inline void MLPutRule(MLINK link, long number, const std::vector& heads = {}) +{ + MLPutFunction(link, "Rule", 2); + MLPutHeads(link, heads); + MLPutLongInteger(link, number); +} + +template +void MLPutRuleTo(MLINK link, T1 t, const T2& name, const std::vector& heads = {}) { MLPutRule(link, name, heads); MLPut(link, t); diff --git a/src/numerics.cpp b/src/numerics.cpp index 7029516b0..249e8f903 100644 --- a/src/numerics.cpp +++ b/src/numerics.cpp @@ -313,21 +313,20 @@ double b22(double p, double m1, double m2, double q) // m1 == m2 with good accuracy if (is_close(m1, m2, EPSTOL)) { answer = -m12 * log(sqr(m1 / q)) * 0.5 + m12 * 0.5; - } - else + } else { /// This zero p limit is good if (fabs(m1) > EPSTOL && fabs(m2) > EPSTOL) { answer = 0.375 * (m12 + m22) - 0.25 * (sqr(m22) * log(sqr(m2 / q)) - sqr(m12) * log(sqr(m1 / q))) / (m22 - m12); - } - else + } else { if (fabs(m1) < EPSTOL) { answer = 0.375 * m22 - 0.25 * m22 * log(sqr(m2 / q)); - } - else { + } else { answer = 0.375 * m12 - 0.25 * m12 * log(sqr(m1 / q)); } + } + } } else { const double b0Save = b0(p, m1, m2, q); diff --git a/templates/LesHouches.in b/templates/LesHouches.in index 177dbb574..946b98630 100644 --- a/templates/LesHouches.in +++ b/templates/LesHouches.in @@ -31,7 +31,7 @@ Block FlexibleSUSYInput 1 125.09 # Mh pole Block SMINPUTS # Standard Model inputs 1 1.279440000e+02 # alpha^(-1) SM MSbar(MZ) - 2 1.166380000e-05 # G_Fermi + 2 1.166378700e-05 # G_Fermi 3 1.184000000e-01 # alpha_s(MZ) SM MSbar 4 9.118760000e+01 # MZ(pole) 5 4.180000000e+00 # mb(mb) SM MSbar diff --git a/templates/edm.cpp.in b/templates/edm.cpp.in new file mode 100644 index 000000000..0053905b8 --- /dev/null +++ b/templates/edm.cpp.in @@ -0,0 +1,630 @@ +// ==================================================================== +// This file is part of FlexibleSUSY. +// +// FlexibleSUSY is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published +// by the Free Software Foundation, either version 3 of the License, +// or (at your option) any later version. +// +// FlexibleSUSY is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with FlexibleSUSY. If not, see +// . +// ==================================================================== + +// File generated at @DateAndTime@ + +/** + * @file @ModelName@_edm.cpp + * + * This file was generated at @DateAndTime@ with FlexibleSUSY + * @FlexibleSUSYVersion@ and SARAH @SARAHVersion@ . + */ + +#include "@ModelName@_edm.hpp" +#include "@ModelName@_mass_eigenstates.hpp" + +#include "numerics2.hpp" +#include "wrappers.hpp" + +#include +#include +#include + +#define INPUTPARAMETER(p) context.model.get_input().p +#define MODELPARAMETER(p) context.model.get_##p() +#define DERIVEDPARAMETER(p) context.model.p() +#define PHASE(p) context.model.get_##p() + +using namespace flexiblesusy; + +namespace { +static constexpr double oneOver16PiSquared = 0.0063325739776461107152; + +/** \brief Helper template for make_array + */ +template +struct make_array_it +{ + using decayed_type = typename std::decay< + typename std::iterator_traits::value_type + >::type; + static constexpr bool use_cast = !std::is_same::value; + + template + static auto iterate( ForwardIterator begin, Args &&...args ) + -> typename std::enable_if< + !use_cast, + std::array + >::type + { + ForwardIterator copy( begin ); + return make_array_it::iterate( ++begin, + std::forward( args )..., *copy ); + } + + template + static auto iterate( ForwardIterator begin, Args &&...args ) + -> typename std::enable_if< + use_cast, + std::array + >::type + { + ForwardIterator copy( begin ); + return make_array_it::iterate( ++begin, + std::forward( args )..., T( *copy ) ); + } +}; + +/** \brief Specialized helper template for make_array */ +template +struct make_array_it +{ + template + static auto iterate( ForwardIterator begin, Args &&...args ) + -> std::array + { + return std::array{ std::forward( args )... }; + } +}; + +/** \class make_array +* \brief Facility for easy construction of \a std::array +* objects. +* \tparam N The size of the to be constructed array +* \tparam T The type of the objects the array should hold +* or void (by default) if the type should be +* inferred. +*/ +template +struct make_array +{ + /** \brief Construct a \a std::array by copying values + * from a range. + * \tparam ForwardIterator The iterator type (usually inferred) + * \param[in] begin The iterator marking the beginning of + * the range to be copied. + * \returns A \a std::array holding the desired objects. + * \warning \a begin must be incrementable sufficiently + * often (\a N times) otherwise the behaviour + * is undefined. + */ + template + static auto iterate( ForwardIterator begin ) + -> std::array::value, + typename std::iterator_traits::value_type, + T + >::type, N> + { + using value_type = typename std::conditional< + std::is_void::value, + typename std::iterator_traits::value_type, + T + >::type; + + return make_array_it::iterate( begin ); + } +}; + +template +auto concatenate( const Array1 &a1, const Array2 &a2 ) +-> std::array< + typename std::common_type< + typename Array1::value_type, + typename Array1::value_type + >::type, + std::tuple_size::value + std::tuple_size::value> +{ + using value_type = typename std::common_type< + typename Array1::value_type, + typename Array1::value_type + >::type; + constexpr auto size = std::tuple_size::value + std::tuple_size::value; + + auto range1 = boost::make_iterator_range( boost::begin(a1), + boost::end(a1) ); + auto range2 = boost::make_iterator_range( boost::begin(a2), + boost::end(a2) ); + auto joined = boost::join( range1, range2 ); + + return make_array::iterate( boost::begin( joined ) ); +} + +/** + * @class IndexBounds + * @brief A class representing multiple (N) index ranges. + * + * N is a non-negative integer. + * The ranges are specified the c++ way; The range begins are inclusive + * and the range ends are exclusive. Misspecified ranges result in + * undefined behaviour! + * The intended use is to iterate over an IndexBounds using the normal + * begin()/end() syntax, which will iterate over every possible combination + * of the indices within their respective ranges. + * The ranges are const! Once they are initialized they cannot be changed. + * Initialization is done like = {{ beg1, beg2, ... }, { end1, end2, ... }} + */ +template struct IndexBounds; + +/** + * @class IndexIterator + * @brief The iterator class used by IndexBounds. + * + * It only fulfils the input iterator requirements! + * The value_type is a const std::array + * containing the current indices. + * Incrementing the iterator results in a new index combination. + */ +template struct IndexIterator +{ + using difference_type = typename std::array::difference_type; + using value_type = unsigned; + using pointer = unsigned *; + using reference = unsigned &; + using iterator_category = std::input_iterator_tag; + + friend class IndexBounds; +private: + const IndexBounds &bounds; + std::array indices; + + IndexIterator(const IndexBounds &b, std::array &&i ) + : bounds(b), indices( std::move(i) ) {} +public: + const std::array &operator*() const { + return indices; + } + + const std::array *operator->() const { + return &indices; + } + + IndexIterator &operator++() { + for (unsigned i = 0; i != N; i++) { + indices[i]++; + if (indices[i] == bounds.indexEnd[i]) { + indices[i] = bounds.indexBegin[i]; + continue; + } + + return *this; + } + + indices = make_array::iterate( bounds.indexEnd ); + return *this; + } + + IndexIterator operator++(int) { + IndexIterator it(*this); + this->operator++(); + return it; + } + + template + friend bool operator==(const IndexIterator &it1, const IndexIterator &it2); +}; + +template +bool operator==(const IndexIterator &it1, const IndexIterator &it2) +{ + return it1.indices == it2.indices; +} + +template +bool operator!=(const IndexIterator &it1, const IndexIterator &it2) +{ + return !(it1 == it2); +} + +template struct IndexBounds { + typedef const std::array indices_type; + + unsigned indexBegin[N]; + unsigned indexEnd[N]; + + typedef IndexIterator const_iterator; + + const_iterator begin() const { + return const_iterator(*this, make_array::iterate( indexBegin )); + } + + const_iterator end() const { + return const_iterator(*this, make_array::iterate( indexEnd )); + } +}; + +template<> struct IndexBounds<0> { + using indices_type = const std::array; + using const_iterator = indices_type *; + + static constexpr indices_type dummyIndex = {}; + + const_iterator begin() const + { return &dummyIndex; } + const_iterator end() const + { return (begin()+1); } +}; +constexpr IndexBounds<0>::indices_type IndexBounds<0>::dummyIndex; + +struct Field {}; + +// Anti fields +template struct anti : public Field +{ + static const unsigned numberOfGenerations = F::numberOfGenerations; + using type = anti; +}; + +template struct anti> { using type = F; }; + +template struct is_anti { static constexpr bool value = false; }; +template struct is_anti> { static constexpr bool value = true; }; + +template struct field_indices; + +template +struct field_indices> +{ + using type = typename field_indices::type; +}; + +@EDM_Fields@ + +/** + * @class EvaluationContext + * @brief Represents an evaluation context. + * + * Actually it simply contains a reference to a model object. + * All computational functions are forwarded to that object, + * e.g. mass calculation functions. + */ +struct EvaluationContext { + @ModelName@_mass_eigenstates& model; ///< The model object. + + /** + * @fn mass + * @brief Returns the mass of a field. + */ + template + double mass( const typename std::enable_if::value, + typename field_indices::type>::type &indices ) const + { return mass::type>( indices ); } + + template + double mass( const typename std::enable_if::value, + typename field_indices::type>::type &indices ) const; + + template + double charge( const typename field_indices::type & ) const; +}; + +/** + * @class DiagramEvaluator + * @brief A template that can calculate disagrams. + * + * This template is completely general. The intended use is + * to have specialisations of the kind: + * DiagramEvaluator + * along with a static member function: + * static double value( EvaluationContext &context ); + * that calculates the contribution to the edm + * for a specific diagram class with specified fields. + */ +template struct DiagramEvaluator; + +double OneLoopFunctionF1C(double); +double OneLoopFunctionF2C(double); +double OneLoopFunctionF1N(double); +double OneLoopFunctionF2N(double); + +@EDM_Diagrams@ + +/** + * @class SingleComponentedVertex + * @brief A vertex whose value can be represented by a single complex number + */ +class SingleComponentedVertex { +private: + std::complex val; ///< The value +public: + SingleComponentedVertex(std::complex v) + : val(v) {} + + /** + * @fn value + * @brief Returns the value of the vertex. + */ + std::complex value() const { + return val; + } + + /** + * @fn isZero + * @brief Tests whether the value is numerically significant + */ + bool isZero() const { + return (is_zero(val.real()) && is_zero(val.imag())); + } +}; + +/** + * @class LeftAndRightComponentedVertex + * @brief A vertex whose value can be represented by two complex numbers + */ +class LeftAndRightComponentedVertex { +private: + std::pair, std::complex> value; ///< The values +public: + LeftAndRightComponentedVertex(const std::complex &left, + const std::complex &right) + : value(left, right) {} + + /** + * @fn left + * @brief Returns the left component of the vertex. + */ + std::complex left() const { + return value.first; + } + + /** + * @fn right + * @brief Returns the left component of the vertex. + */ + std::complex right() const { + return value.second; + } + + /** + * @fn isZero + * @brief Tests whether the values are numerically significant + */ + bool isZero() const { + return (is_zero(value.first.real()) && is_zero(value.first.imag()) && + is_zero(value.second.real()) && is_zero(value.second.imag())); + } +}; + +/** + * @class VertexFunctionData + * @brief Vertex data for a vertex with the fields specified by F... + * + * All elements in F... have to be publicly derived from Field. + * The data includes the type of the vertex, whether F... is in + * canonical order or permuted as well as index bounds referring to + * the specific fields (generation and colour indices, etc). + * The intention is to have a program generate VertexFunctionData<> + * specialisations that are then used by the VertexFunction<> template. + */ +template struct VertexFunctionData; + +/** + * @class VertexFunction + * @brief A template that represents a vertex with open field indices. + * + * All elements in F... have to be publicly derived from Field. + * To obtain a conrete value, use the static member function vertex() + * along with the desired field indices. + */ +template class VertexFunction { + using Data = VertexFunctionData; + using bounds_type = typename std::decay::type; +public: + using vertex_type = typename Data::vertex_type; + using indices_type = typename decltype(Data::index_bounds)::indices_type; + + static constexpr bounds_type index_bounds = Data::index_bounds; + + /** + * @fn fieldIndices + * @brief Returns the subset of indices belonging to + * a field at a given index. + * @param indices The complete set of indices for all fields + * + * The field indexing is in the same order as the template arguments + * and starts at 0. + */ + template + static + std::array + fieldIndices(const indices_type &indices) + { + constexpr auto length = Data::fieldIndexStart[fieldIndex+1] - Data::fieldIndexStart[fieldIndex]; + constexpr auto offset = Data::fieldIndexStart[fieldIndex]; + auto begin = indices.begin() + offset; + return make_array::iterate(begin); + } + + /** + * @fn vertex + * @brief Calculates the vertex for a given set of field indices. + * @param indices The field indices + * @param context The evaluation context + */ + static vertex_type vertex(const indices_type &indices, const EvaluationContext &context); +}; + +@EDM_VertexFunctionData@ +} + +@EDM_InterfaceFunctionDefinitions@ + +namespace { +/** +* @defgroup LoopFunctions +* @brief The loop functions necessary for edm one-loop calculations. +* +* These are OneLoopFunctionA(), OneLoopFunctionB() +* as specified in arXiv:0808.1819 +*/ + +double OneLoopFunctionA(double r) +{ + if( is_zero(r) ) + return std::numeric_limits::infinity(); + + const double d = r - 1.0; + + if( std::abs(d) < 0.15 ) { + return (-0.33333333333333333333 + + 0.25000000000000000000 * d - + 0.20000000000000000000 * d * d + + 0.16666666666666666667 * d * d * d - + 0.14285714285714285714 * d * d * d * d + + 0.12500000000000000000 * d * d * d * d * d - + 0.11111111111111111111 * d * d * d * d * d * d + + 0.10000000000000000000 * d * d * d * d * d * d * d - + 0.090909090909090909091 * d * d * d * d * d * d * d * d + + 0.083333333333333333333 * d * d * d * d * d * d * d * d * d - + 0.076923076923076923077 * d * d * d * d * d * d * d * d * d * d ); + } + + return 1.0 / (2.0 * d * d) * (3.0 - r - 2.0 * std::log(r) / d); +} + +double OneLoopFunctionB(double r) +{ + if( is_zero(r) ) + return 1.0/2.0; + + const double d = r - 1.0; + + if( std::abs(d) < 0.15 ) + return (0.16666666666666666667 - + 0.083333333333333333333 * d + + 0.050000000000000000000 * d * d - + 0.033333333333333333333 * d * d * d + + 0.023809523809523809524 * d * d * d * d - + 0.017857142857142857143 * d * d * d * d * d + + 0.013888888888888888889 * d * d * d * d * d * d - + 0.011111111111111111111 * d * d * d * d * d * d * d + + 0.0090909090909090909091 * d * d * d * d * d * d * d * d - + 0.0075757575757575757576 * d * d * d * d * d * d * d * d * d + + 0.0064102564102564102564 * d * d * d * d * d * d * d * d * d * d ); + + return 1.0 / (2.0 * d * d) * (1.0 + r - 2.0 * r * std::log(r) / d); +} + +@EDM_Definitions@ + +// PhotonEmitter is the fermion +template +double DiagramEvaluator, + EDMField, PhotonEmitter, ExchangeField +>::value(const typename field_indices::type &indices, EvaluationContext &context) +{ + double res = 0.0; + + using FermionVertex = VertexFunction< + typename anti::type, + PhotonEmitter, + ExchangeField + >; + + constexpr auto indexBounds = FermionVertex::index_bounds; + for( auto indexIt = indexBounds.begin(); indexIt != indexBounds.end(); ++indexIt ) + { + auto edmFieldIndices = FermionVertex::template fieldIndices<0>(*indexIt); + + if( edmFieldIndices != indices ) + continue; + + auto photonEmitterIndices = FermionVertex::template fieldIndices<1>(*indexIt); + auto exchangeFieldIndices = FermionVertex::template fieldIndices<2>(*indexIt); + auto vertex = FermionVertex::vertex(*indexIt, context); + + auto photonEmitterMass = context.mass( photonEmitterIndices ); + auto exchangeFieldMass = context.mass( exchangeFieldIndices ); + + double photonEmitterCharge = context.charge( photonEmitterIndices ); + + constexpr double numericFactor = oneOver16PiSquared; + double massFactor = photonEmitterMass/(exchangeFieldMass * exchangeFieldMass); + double couplingFactor = (std::conj( vertex.right() ) * vertex.left()).imag(); + + double massRatioSquared = photonEmitterMass/exchangeFieldMass; + massRatioSquared *= massRatioSquared; + + double loopFactor = photonEmitterCharge * OneLoopFunctionA( massRatioSquared ); + + double contribution = numericFactor * massFactor * couplingFactor * loopFactor; + + res += contribution; + } + + return res; +} + +// ExchangeField is fermion +template +double DiagramEvaluator, + EDMField, PhotonEmitter, ExchangeField +>::value(const typename field_indices::type &indices, EvaluationContext &context) +{ + double res = 0.0; + + using FermionVertex = VertexFunction< + typename anti::type, + PhotonEmitter, + ExchangeField + >; + + constexpr auto indexBounds = FermionVertex::index_bounds; + for( auto indexIt = indexBounds.begin(); indexIt != indexBounds.end(); ++indexIt ) + { + auto edmFieldIndices = FermionVertex::template fieldIndices<0>(*indexIt); + + if( edmFieldIndices != indices ) + continue; + + auto photonEmitterIndices = FermionVertex::template fieldIndices<1>(*indexIt); + auto exchangeFieldIndices = FermionVertex::template fieldIndices<2>(*indexIt); + auto vertex = FermionVertex::vertex(*indexIt, context); + + auto photonEmitterMass = context.mass( photonEmitterIndices ); + auto exchangeFieldMass = context.mass( exchangeFieldIndices ); + + double photonEmitterCharge = context.charge( photonEmitterIndices ); + + constexpr double numericFactor = oneOver16PiSquared; + double massFactor = exchangeFieldMass/(photonEmitterMass * photonEmitterMass); + double couplingFactor = (std::conj( vertex.right() ) * vertex.left()).imag(); + + double massRatioSquared = exchangeFieldMass / photonEmitterMass; + massRatioSquared *= massRatioSquared; + + double loopFactor = photonEmitterCharge * OneLoopFunctionB( massRatioSquared ); + + double contribution = numericFactor * massFactor * couplingFactor * loopFactor; + + res += contribution; + } + + return res; +} +} diff --git a/templates/edm.hpp.in b/templates/edm.hpp.in new file mode 100644 index 000000000..65712b80f --- /dev/null +++ b/templates/edm.hpp.in @@ -0,0 +1,37 @@ +// ==================================================================== +// This file is part of FlexibleSUSY. +// +// FlexibleSUSY is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published +// by the Free Software Foundation, either version 3 of the License, +// or (at your option) any later version. +// +// FlexibleSUSY is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with FlexibleSUSY. If not, see +// . +// ==================================================================== + +// File generated at @DateAndTime@ + +/** + * @file @ModelName@_edm.hpp + * + * This file was generated at @DateAndTime@ with FlexibleSUSY + * @FlexibleSUSYVersion@ and SARAH @SARAHVersion@ . + */ + +#ifndef @ModelName@_EDM_H +#define @ModelName@_EDM_H + +namespace flexiblesusy { +class @ModelName@_mass_eigenstates; + +@EDM_InterfaceFunctionPrototypes@ +} // namespace flexiblesusy + +#endif diff --git a/templates/librarylink.m.in b/templates/librarylink.m.in index bcb236e3c..8f40df0f9 100644 --- a/templates/librarylink.m.in +++ b/templates/librarylink.m.in @@ -64,7 +64,7 @@ fsDefaultSettings = { fsDefaultSMParameters = { alphaEmMZ -> 1/127.916, (* SMINPUTS[1] *) - GF -> 1.16637*^-5, (* SMINPUTS[2] *) + GF -> 1.1663787*^-5, (* SMINPUTS[2] *) alphaSMZ -> 0.1184, (* SMINPUTS[3] *) MZ -> 91.1876, (* SMINPUTS[4] *) mbmb -> 4.18, (* SMINPUTS[5] *) diff --git a/templates/module.mk b/templates/module.mk index 86f042eb7..8dd9efeea 100644 --- a/templates/module.mk +++ b/templates/module.mk @@ -4,6 +4,8 @@ MODNAME := templates BASE_TEMPLATES := \ $(DIR)/a_muon.hpp.in \ $(DIR)/a_muon.cpp.in \ + $(DIR)/edm.hpp.in \ + $(DIR)/edm.cpp.in \ $(DIR)/convergence_tester.hpp.in \ $(DIR)/ewsb_solver.hpp.in \ $(DIR)/ewsb_solver_interface.hpp.in \ diff --git a/templates/module.mk.in b/templates/module.mk.in index 0c82caa7f..2e4703322 100644 --- a/templates/module.mk.in +++ b/templates/module.mk.in @@ -34,6 +34,7 @@ WITH_$(MODNAME) := yes LIB@CLASSNAME@_SRC := \ $(DIR)/@CLASSNAME@_a_muon.cpp \ + $(DIR)/@CLASSNAME@_edm.cpp \ $(DIR)/@CLASSNAME@_effective_couplings.cpp \ $(DIR)/@CLASSNAME@_info.cpp \ $(DIR)/@CLASSNAME@_input_parameters.cpp \ @@ -62,6 +63,7 @@ LL@CLASSNAME@_MMA := \ LIB@CLASSNAME@_HDR := \ $(DIR)/@CLASSNAME@_a_muon.hpp \ $(DIR)/@CLASSNAME@_convergence_tester.hpp \ + $(DIR)/@CLASSNAME@_edm.hpp \ $(DIR)/@CLASSNAME@_effective_couplings.hpp \ $(DIR)/@CLASSNAME@_ewsb_solver.hpp \ $(DIR)/@CLASSNAME@_ewsb_solver_interface.hpp \ diff --git a/templates/observables.cpp.in b/templates/observables.cpp.in index c0dc44365..889fbea68 100644 --- a/templates/observables.cpp.in +++ b/templates/observables.cpp.in @@ -21,6 +21,7 @@ #include "@ModelName@_observables.hpp" #include "@ModelName@_mass_eigenstates.hpp" #include "@ModelName@_a_muon.hpp" +#include "@ModelName@_edm.hpp" #include "@ModelName@_effective_couplings.hpp" #include "gm2calc_interface.hpp" #include "eigen_utils.hpp" diff --git a/test/module.mk b/test/module.mk index 33277bb4c..b15186f9f 100644 --- a/test/module.mk +++ b/test/module.mk @@ -55,6 +55,31 @@ TEST_SH := \ $(DIR)/test_run_all_spectrum_generators.sh \ $(DIR)/test_space_dir.sh +TEST_META := \ + $(DIR)/test_BetaFunction.m \ + $(DIR)/test_CConversion.m \ + $(DIR)/test_Constraint.m \ + $(DIR)/test_EWSB.m \ + $(DIR)/test_LoopFunctions.m \ + $(DIR)/test_MSSM_2L_analytic.m \ + $(DIR)/test_MSSM_2L_mt.m \ + $(DIR)/test_MSSM_2L_yb_softsusy.m \ + $(DIR)/test_MSSM_2L_yt.m \ + $(DIR)/test_MSSM_2L_yt_loopfunction.m \ + $(DIR)/test_MSSM_2L_yt_softsusy.m \ + $(DIR)/test_Parameters.m \ + $(DIR)/test_ReadSLHA.m \ + $(DIR)/test_RGIntegrator.m \ + $(DIR)/test_SelfEnergies.m \ + $(DIR)/test_SemiAnalytic.m \ + $(DIR)/test_TextFormatting.m \ + $(DIR)/test_THDM_threshold_corrections.m \ + $(DIR)/test_THDM_threshold_corrections_gauge.m \ + $(DIR)/test_ThreeLoopQCD.m \ + $(DIR)/test_ThresholdCorrections.m \ + $(DIR)/test_TreeMasses.m \ + $(DIR)/test_Vertices.m + ifneq ($(findstring lattice,$(SOLVERS)),) TEST_SRC += endif # ifneq($(findstring lattice,$(SOLVERS)),) @@ -398,6 +423,8 @@ endif ifeq ($(WITH_CMSSMCPV),yes) TEST_SRC += \ $(DIR)/test_CMSSMCPV_ewsb.cpp +TEST_META += \ + $(DIR)/test_CMSSMCPV_librarylink.m endif ifeq ($(WITH_CMSSMCPV) $(WITH_CMSSM),yes yes) @@ -504,31 +531,6 @@ TEST_META += \ $(DIR)/test_HGTHDM_THDM_threshold_corrections_scale_invariance.m endif -TEST_META := \ - $(DIR)/test_BetaFunction.m \ - $(DIR)/test_CConversion.m \ - $(DIR)/test_Constraint.m \ - $(DIR)/test_EWSB.m \ - $(DIR)/test_LoopFunctions.m \ - $(DIR)/test_MSSM_2L_analytic.m \ - $(DIR)/test_MSSM_2L_mt.m \ - $(DIR)/test_MSSM_2L_yb_softsusy.m \ - $(DIR)/test_MSSM_2L_yt.m \ - $(DIR)/test_MSSM_2L_yt_loopfunction.m \ - $(DIR)/test_MSSM_2L_yt_softsusy.m \ - $(DIR)/test_Parameters.m \ - $(DIR)/test_ReadSLHA.m \ - $(DIR)/test_RGIntegrator.m \ - $(DIR)/test_SelfEnergies.m \ - $(DIR)/test_SemiAnalytic.m \ - $(DIR)/test_TextFormatting.m \ - $(DIR)/test_THDM_threshold_corrections.m \ - $(DIR)/test_THDM_threshold_corrections_gauge.m \ - $(DIR)/test_ThreeLoopQCD.m \ - $(DIR)/test_ThresholdCorrections.m \ - $(DIR)/test_TreeMasses.m \ - $(DIR)/test_Vertices.m - ifeq ($(WITH_SM),yes) TEST_META += \ $(DIR)/test_SM_3loop_beta.m diff --git a/test/test_CMSSMCPV_librarylink.in.spc b/test/test_CMSSMCPV_librarylink.in.spc new file mode 100644 index 000000000..93cee4b1a --- /dev/null +++ b/test/test_CMSSMCPV_librarylink.in.spc @@ -0,0 +1,54 @@ +Block MODSEL # Select model +# 12 1000 # DRbar parameter output scale (GeV) +Block FlexibleSUSY + 0 1.000000000e-04 # precision goal + 1 0 # max. iterations (0 = automatic) + 2 0 # algorithm (0 = two_scale, 1 = lattice) + 3 0 # calculate SM pole masses + 4 2 # pole mass loop order + 5 2 # EWSB loop order + 6 3 # beta-functions loop order + 7 2 # threshold corrections loop order + 8 1 # Higgs 2-loop corrections O(alpha_t alpha_s) + 9 1 # Higgs 2-loop corrections O(alpha_b alpha_s) + 10 1 # Higgs 2-loop corrections O((alpha_t + alpha_b)^2) + 11 1 # Higgs 2-loop corrections O(alpha_tau^2) + 12 0 # force output + 13 1 # Top quark 2-loop corrections QCD + 14 1.000000000e-11 # beta-function zero threshold + 15 1 # calculate observables (a_muon, ...) + 16 0 # force positive majorana masses + 17 0 # pole mass renormalization scale (0 = SUSY scale) +Block FlexibleSUSYInput + 0 0.00729735 # alpha_em(0) + 1 125.09 # Mh pole +Block SMINPUTS # Standard Model inputs + 1 1.279340000e+02 # alpha^(-1) SM MSbar(MZ) + 2 1.166370000e-05 # G_Fermi + 3 1.176000000e-01 # alpha_s(MZ) SM MSbar + 4 9.118760000e+01 # MZ(pole) + 5 4.200000000e+00 # mb(mb) SM MSbar + 6 1.733000000e+02 # mtop(pole) + 7 1.777000000e+00 # mtau(pole) + 8 0.000000000e+00 # mnu3(pole) + 9 80.404 # MW pole + 11 5.109989020e-04 # melectron(pole) + 12 0.000000000e+00 # mnu1(pole) + 13 1.056583570e-01 # mmuon(pole) + 14 0.000000000e+00 # mnu2(pole) + 21 4.750000000e-03 # md(2 GeV) MS-bar + 22 2.400000000e-03 # mu(2 GeV) MS-bar + 23 1.040000000e-01 # ms(2 GeV) MS-bar + 24 1.270000000e+00 # mc(mc) MS-bar +Block MINPAR # Input parameters + 1 125. # m0 + 2 500. # m12 + 3 10. # TanBeta + 4 1 # SignMu + 5 0 # Azero +Block IMMINPAR + 2 10 # Imm12 + 4 0 # SinPhiMu + 5 10 # ImAzero +Block EXTPAR + 100 0.1 # etaInput diff --git a/test/test_CMSSMCPV_librarylink.m b/test/test_CMSSMCPV_librarylink.m new file mode 100644 index 000000000..79e2bcc77 --- /dev/null +++ b/test/test_CMSSMCPV_librarylink.m @@ -0,0 +1,153 @@ +Needs["TestSuite`", "TestSuite.m"]; +Needs["ReadSLHA`", "ReadSLHA.m"]; + +Get["models/CMSSMCPV/CMSSMCPV_librarylink.m"]; + +Print["Comparing SLHA output to LibraryLink output"]; + +settings = { + precisionGoal -> 0.0001, + maxIterations -> 0, + calculateStandardModelMasses -> 0, + poleMassLoopOrder -> 2, + ewsbLoopOrder -> 2, + betaFunctionLoopOrder -> 3, + thresholdCorrectionsLoopOrder -> 2, + higgs2loopCorrectionAtAs -> 1, + higgs2loopCorrectionAbAs -> 1, + higgs2loopCorrectionAtAt -> 1, + higgs2loopCorrectionAtauAtau -> 1, + forceOutput -> 0, + topPoleQCDCorrections -> 1, + betaZeroThreshold -> 1.*10^-11, + forcePositiveMasses -> 0, + poleMassScale -> 0., + parameterOutputScale -> 0. +}; + +smInputs = { + alphaEmMZ -> 1./127.934, + GF -> 0.0000116637, + alphaSMZ -> 0.1176, + MZ -> 91.1876, + mbmb -> 4.2, + Mt -> 173.3, + Mtau -> 1.777, + Mv3 -> 0, + MW -> 80.404, + Me -> 0.000510998902, + Mv1 -> 0., + Mm -> 0.105658357, + Mv2 -> 0., + md2GeV -> 0.00475, + mu2GeV -> 0.0024, + ms2GeV -> 0.104, + mcmc -> 1.27 +}; + +handle = FSCMSSMCPVOpenHandle[ + fsSettings -> settings, + fsSMParameters -> smInputs, + fsModelParameters -> { + m0 -> 125, m12 -> 500, TanBeta -> 10, CosPhiMu -> 1, Azero -> 0, + Imm12 -> 10, SinPhiMu -> 0, ImAzero -> 10, etaInput -> 0.1 + } +]; + +specML = CMSSMCPV /. FSCMSSMCPVCalculateSpectrum[handle]; +obsML = FSCMSSMCPVCalculateObservables[handle]; +probML = FSCMSSMCPVGetProblems[handle]; +warnML = FSCMSSMCPVGetWarnings[handle]; + +FSCMSSMCPVCloseHandle[handle]; + +inputFile = "test/test_CMSSMCPV_librarylink.in.spc"; +outputFile = "test/test_CMSSMCPV_librarylink.out.spc"; +cmd = "models/CMSSMCPV/run_CMSSMCPV.x --slha-input-file=" <> inputFile <> + " --slha-output-file=" <> outputFile; + +Run["rm -f " <> outputFile]; +Run[cmd]; +slhaStr = Import[outputFile, "String"]; + +parameters = { + {g1, {0} , {GAUGE, 1}}, + {g2, {0} , {GAUGE, 2}}, + {g3, {0} , {GAUGE, 3}}, + {Yu, {3, 3}, Yu}, + {Yd, {3, 3}, Yd}, + {Ye, {3, 3}, Ye}, + {Tu, {3, 3}, Tu}, + {Td, {3, 3}, Td}, + {Te, {3, 3}, Te}, + {ImTu, {3, 3}, ImTu}, + {ImTd, {3, 3}, ImTd}, + {ImTe, {3, 3}, ImTe}, + {mq2, {3, 3}, MSQ2}, + {ml2, {3, 3}, MSL2}, + {mu2, {3, 3}, MSU2}, + {md2, {3, 3}, MSD2}, + {me2, {3, 3}, MSE2}, + {Immq2, {3, 3}, ImMSQ2}, + {Imml2, {3, 3}, ImMSL2}, + {Immu2, {3, 3}, ImMSU2}, + {Immd2, {3, 3}, ImMSD2}, + {Imme2, {3, 3}, ImMSE2}, + {mHu2, {0}, {MSOFT, 22}}, + {mHd2, {0}, {MSOFT, 21}}, + {vu, {0}, {HMIX, 103}}, + {vd, {0}, {HMIX, 102}}, + {Mu, {0}, {HMIX, 1}}, + {BMu, {0}, {HMIX, 101}}, + {ImMu, {0}, {ImHMIX, 1}}, + {ImBMu, {0}, {ImHMIX, 101}}, + {CpHPP1, {0}, {EFFHIGGSCOUPLINGS, 25, 22, 22}}, + {CpHPP2, {0}, {EFFHIGGSCOUPLINGS, 35, 22, 22}}, + {CpHPP3, {0}, {EFFHIGGSCOUPLINGS, 36, 22, 22}}, + {CpHGG1, {0}, {EFFHIGGSCOUPLINGS, 25, 21, 21}}, + {CpHGG2, {0}, {EFFHIGGSCOUPLINGS, 35, 21, 21}}, + {CpHGG3, {0}, {EFFHIGGSCOUPLINGS, 36, 21, 21}}, + {aMuon, {0}, {FlexibleSUSYLowEnergy, 21}}, + {EDM[Fe[0]], {0}, {FlexibleSUSYLowEnergy, 23}}, + {EDM[Fe[1]], {0}, {FlexibleSUSYLowEnergy, 24}}, + {EDM[Fe[2]], {0}, {FlexibleSUSYLowEnergy, 25}} +}; + +slhaData = ReadSLHAString[slhaStr, parameters]; + +delta = 1*^-4; + +TestCloseRel[g1 * Sqrt[5/3] /. slhaData, g1 /. specML, delta]; +TestCloseRel[g2 /. slhaData, g2 /. specML, delta]; +TestCloseRel[g3 /. slhaData, g3 /. specML, delta]; +TestCloseRel[Abs[Yu] /. slhaData, Abs[Yu] /. specML, delta]; +TestCloseRel[Abs[Yd] /. slhaData, Abs[Yd] /. specML, delta]; +TestCloseRel[Abs[Ye] /. slhaData, Abs[Ye] /. specML, delta]; +TestCloseRel[Abs[Tu + I ImTu] /. slhaData, Abs[T[Yu]] /. specML, delta]; +TestCloseRel[Abs[Td + I ImTd] /. slhaData, Abs[T[Yd]] /. specML, delta]; +TestCloseRel[Abs[Te + I ImTe] /. slhaData, Abs[T[Ye]] /. specML, delta]; +TestCloseRel[Abs[mq2 + I Immq2] /. slhaData, Abs[mq2] /. specML, delta]; +TestCloseRel[Abs[ml2 + I Imml2] /. slhaData, Abs[ml2] /. specML, delta]; +TestCloseRel[Abs[mu2 + I Immu2] /. slhaData, Abs[mu2] /. specML, delta]; +TestCloseRel[Abs[md2 + I Immd2] /. slhaData, Abs[md2] /. specML, delta]; +TestCloseRel[Abs[me2 + I Imme2] /. slhaData, Abs[me2] /. specML, delta]; +TestCloseRel[mHu2 /. slhaData, mHu2 /. specML, delta]; +TestCloseRel[mHd2 /. slhaData, mHd2 /. specML, delta]; +TestCloseRel[vu /. slhaData, vu /. specML, delta]; +TestCloseRel[vd /. slhaData, vd /. specML, delta]; +TestCloseRel[Mu + I ImMu /. slhaData, \[Mu] /. specML, delta]; +TestCloseRel[BMu + I ImBMu /. slhaData, B[\[Mu]] /. specML, delta]; + +delta = 1*^-5; + +TestCloseRel[{CpHPP1, CpHPP2, CpHPP3} /. slhaData, Abs[FlexibleSUSYObservable`CpHiggsPhotonPhoton /. obsML], delta]; +TestCloseRel[{CpHGG1, CpHGG2, CpHGG3} /. slhaData, Abs[FlexibleSUSYObservable`CpHiggsGluonGluon] /. obsML, delta]; +TestCloseRel[aMuon /. slhaData, FlexibleSUSYObservable`aMuon /. obsML, delta]; +TestCloseRel[EDM[Fe[0]] /. slhaData, FlexibleSUSYObservable`EDM[Fe[0]] /. obsML, delta]; +TestCloseRel[EDM[Fe[1]] /. slhaData, FlexibleSUSYObservable`EDM[Fe[1]] /. obsML, delta]; +TestCloseRel[EDM[Fe[2]] /. slhaData, FlexibleSUSYObservable`EDM[Fe[2]] /. obsML, delta]; + +TestEquality[probML, {}]; +TestEquality[warnML, {}]; + +PrintTestSummary[]; diff --git a/test/test_CMSSM_librarylink.m b/test/test_CMSSM_librarylink.m index 125bf9d50..7f08c3814 100644 --- a/test/test_CMSSM_librarylink.m +++ b/test/test_CMSSM_librarylink.m @@ -94,7 +94,8 @@ {CpHGG1, {0}, {EFFHIGGSCOUPLINGS, 25, 21, 21}}, {CpHGG2, {0}, {EFFHIGGSCOUPLINGS, 35, 21, 21}}, {CpAPP, {0}, {EFFHIGGSCOUPLINGS, 36, 22, 22}}, - {CpAGG, {0}, {EFFHIGGSCOUPLINGS, 36, 21, 21}} + {CpAGG, {0}, {EFFHIGGSCOUPLINGS, 36, 21, 21}}, + {aMuon, {0}, {FlexibleSUSYLowEnergy, 21}} }; slhaData = ReadSLHAString[slhaStr, parameters]; @@ -128,6 +129,7 @@ TestCloseRel[{CpHGG1, CpHGG2} /. slhaData, Abs[FlexibleSUSYObservable`CpHiggsGluonGluon] /. obsML, delta]; TestCloseRel[CpAPP /. slhaData, Abs[FlexibleSUSYObservable`CpPseudoScalarPhotonPhoton /. obsML], delta]; TestCloseRel[CpAGG /. slhaData, Abs[FlexibleSUSYObservable`CpPseudoScalarGluonGluon /. obsML], delta]; +TestCloseRel[aMuon /. slhaData, FlexibleSUSYObservable`aMuon /. obsML, delta]; TestEquality[probML, {}]; TestEquality[warnML, {}];