Skip to content

Commit

Permalink
And... We compile\!
Browse files Browse the repository at this point in the history
  • Loading branch information
@iolojz
iolojz committed Sep 3, 2017
1 parent 55bca37 commit 7b16e4a
Show file tree
Hide file tree
Showing 5 changed files with 72 additions and 44 deletions.
16 changes: 13 additions & 3 deletions meta/AMuon.m
View file
Expand Up @@ -52,7 +52,17 @@

Return[False];
]


CreateCalculateMuonPoleMass[] := "model.calculate_M" <>
CXXDiagrams`CXXNameOfField[GetMuon[]] <>
"_pole();"
CreateMuonPhysicalMass[] := "return context.model.get_physical().M" <>
CXXDiagrams`CXXNameOfField[GetMuon[]] <>
"(" <> If[GetMuonIndex[] =!= Null,
" " <> ToString @ GetMuonIndex[] <> " ",
""] <>
");";

CreateCalculation[gTaggedDiagrams_List] :=
Module[{muon = GetMuon[], muonIndex = GetMuonIndex[],
calculation,numberOfIndices},
Expand Down Expand Up @@ -132,7 +142,7 @@

GetQED2L[] :=
Module[{muonIndex = GetMuonIndex[],
numberOfIndices = CXXDiagrams`NumberOfFieldIndices[GetMuon]},
numberOfIndices = CXXDiagrams`NumberOfFieldIndices[GetMuon[]]},
"const field_indices<Muon>::type muonIndices = {" <>
If[muonIndex =!= Null,
" " <> ToString @ muonIndex <>
Expand All @@ -148,7 +158,7 @@

"const double MSUSY = Abs(get_MSUSY(context.model));\n" <>
"const double m_muon = muonPhysicalMass(context);\n" <>
"const double alpha_em = Sqr(context.charge<Muon>( muonIndices ))/(4*Pi);\n" <>
"const double alpha_em = Sqr(context.electric_charge<Muon>( muonIndices ))/(4*Pi);\n" <>
"const double qed_2L = alpha_em/(4*Pi) * 16 * FiniteLog(m_muon/MSUSY);\n\n" <>
"return qed_2L;"
]
Expand Down
7 changes: 2 additions & 5 deletions meta/CXXDiagrams.m
View file
Expand Up @@ -44,7 +44,7 @@
"conj"];
LorentzConjugate[field_] := SARAH`AntiField[field]

CreateFields[specialFields_List] :=
CreateFields[] :=
Module[{fields},
fields = TreeMasses`GetParticles[];

Expand All @@ -58,11 +58,8 @@
CXXNameOfField[LorentzConjugate[#]] <> ";\n"] <>
"};\n" &) /@ fields, "\n"] <> "\n\n" <>

"// Special fields\n" <>
"// Named fields\n" <>
"using Photon = " <> CXXNameOfField[SARAH`Photon] <> ";\n\n" <>
StrinJoin @ Riffle[
("using " <> #[[1]] <> " = " <> CXXNameOfField @ #[[2]] <> ";" &) /@
Cases[specialFields,Except[{"Photon",_}]],"\n"] <> "\n\n" <>

"// Fields that are their own Lorentz conjugates.\n" <>
StringJoin @ Riffle[
Expand Down
27 changes: 17 additions & 10 deletions meta/FlexibleSUSY.m
View file
Expand Up @@ -1822,11 +1822,11 @@ corresponding tadpole is real or imaginary (only in models with CP
];

(* Write the CXXDiagrams c++ files *)
WriteCXXDiagramClass[specialFields_List,vertices_List,massMatrices_,files_List] :=
WriteCXXDiagramClass[vertices_List,massMatrices_,files_List] :=
Module[{fields, nPointFunctions, vertexRules, vertexData, cxxVertices, massFunctions},
vertexRules = CXXDiagrams`VertexRulesForVertices[vertices,massMatrices];

fields = CXXDiagrams`CreateFields[specialFields];
fields = CXXDiagrams`CreateFields[];
vertexData = StringJoin @ Riffle[CXXDiagrams`CreateVertexData[#,vertexRules] &
/@ DeleteDuplicates[vertices],
"\n\n"];
Expand Down Expand Up @@ -1870,24 +1870,31 @@ corresponding tadpole is real or imaginary (only in models with CP
(* Write the AMuon c++ files *)
WriteAMuonClass[vertexRules_List, files_List] :=
Module[{graphs,diagrams,vertices,
muonPoleMass,
muonPhysicalMass,
calculation,
getMSUSY, getQED2L},
graphs = AMuon`ContributingGraphs[];
diagrams = AMuon`ContributingDiagramsForGraph /@ graphs;

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

muonPoleMass = AMuon`CreateCalculateMuonPoleMass[];
muonPhysicalMass = AMuon`CreateMuonPhysicalMass[];
calculation = AMuon`CreateCalculation @ Transpose[{graphs,diagrams}];

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

WriteOut`ReplaceInFiles[files,
{"@AMuon_Calculation@" -> TextFormatting`IndentText[calculation],
"@AMuon_GetMSUSY@" -> IndentText[WrapLines[getMSUSY]],
"@AMuon_QED_2L@" -> IndentText[WrapLines[getQED2L]],
Sequence @@ GeneralReplacementRules[]
}];
{"@AMuon_MuonField@" -> CXXDiagrams`CXXNameOfField[AMuon`GetMuon[]],
"@AMuon_CalculateMuonPoleMass@" -> muonPoleMass,
"@AMuon_MuonPhysicalMass@" -> TextFormatting`IndentText[muonPhysicalMass],
"@AMuon_Calculation@" -> TextFormatting`IndentText[calculation],
"@AMuon_GetMSUSY@" -> IndentText[WrapLines[getMSUSY]],
"@AMuon_QED_2L@" -> IndentText[WrapLines[getQED2L]],
Sequence @@ GeneralReplacementRules[]
}];

vertices
];
Expand Down Expand Up @@ -2968,7 +2975,7 @@ corresponding tadpole is real or imaginary (only in models with CP

MakeFlexibleSUSY[OptionsPattern[]] :=
Module[{nPointFunctions, runInputFile, initialGuesserInputFile,
edm2Vertices,aMuonVertices,edm2NPointFunctions,cxxVertexRules,
edm2Vertices,aMuonVertices,
gmm2Vertices = {}, edmFields,
susyBetaFunctions, susyBreakingBetaFunctions,
numberOfSusyParameters, anomDim,
Expand Down Expand Up @@ -3847,8 +3854,8 @@ corresponding tadpole is real or imaginary (only in models with CP
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"}]}}];
{{FileNameJoin[{$flexiblesusyTemplateDir, "cxx_diagrams.hpp.in"}],
FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> "_cxx_diagrams.hpp"}]}}];



Expand Down
58 changes: 36 additions & 22 deletions templates/a_muon2.cpp.in
View file
Expand Up @@ -19,7 +19,7 @@
// File generated at @DateAndTime@

/**
* @file @ModelName@_a_muon.cpp
* @file @ModelName@_a_muon2.cpp
*
* This file was generated at @DateAndTime@ with FlexibleSUSY
* @FlexibleSUSYVersion@ and SARAH @SARAHVersion@ .
Expand All @@ -38,6 +38,8 @@
using namespace flexiblesusy;
using namespace cxx_diagrams;

using Muon = @AMuon_MuonField@;

namespace {
static constexpr double oneOver16PiSquared = 0.0063325739776461107152;

Expand All @@ -61,7 +63,7 @@ template<class PhotonEmitter, class ExchangeParticle>
struct AMuonVertexCorrectionSF
{
static double value(const typename field_indices<Muon>::type &indices,
EvaluationContext &context);
const EvaluationContext &context);
};

/**
Expand All @@ -84,11 +86,9 @@ template<class PhotonEmitter, class ExchangeParticle>
struct AMuonVertexCorrectionFS
{
static double value(const typename field_indices<Muon>::type &indices,
EvaluationContext &context);
const EvaluationContext &context);
};
}

namespace {
/**
* @defgroup LoopFunctions LoopFunctions
* @brief The loop functions necessary for \f$a_\mu\f$ one-loop calculations.
Expand Down Expand Up @@ -300,9 +300,13 @@ std::pair<double,double> vary_scale(const @ModelName@_mass_eigenstates& model)
return std::make_pair(*(minmax.first), *(minmax.second));
}

double muonPhysicalMass( const EvaluationContext &context )
{
@AMuon_MuonPhysicalMass@
}
} // anonymous namespace

double @ModelName@_a_muon::calculate_a_muon(const @ModelName@_mass_eigenstates& model_)
double @ModelName@_a_muon2::calculate_a_muon(const @ModelName@_mass_eigenstates& model_)
{
@ModelName@_mass_eigenstates model(model_);

Expand All @@ -316,10 +320,16 @@ double @ModelName@_a_muon::calculate_a_muon(const @ModelName@_mass_eigenstates&
return std::numeric_limits<double>::quiet_NaN();
}

double m_muon_pole = muonPhysicalMass( EvaluationContext{ model } );
if (m_muon_pole == 0.0) {
model.solve_ewsb();
@AMuon_CalculateMuonPoleMass@
}

return calculate_a_muon_impl(model);
}

double @ModelName@_a_muon::calculate_a_muon_uncertainty(const @ModelName@_mass_eigenstates& model_)
double @ModelName@_a_muon2::calculate_a_muon_uncertainty(const @ModelName@_mass_eigenstates& model_)
{
@ModelName@_mass_eigenstates model(model_);

Expand All @@ -339,6 +349,7 @@ double @ModelName@_a_muon::calculate_a_muon_uncertainty(const @ModelName@_mass_e
return delta_amu_scale;
}

namespace {
double get_QED_2L(EvaluationContext& context)
{
@AMuon_QED_2L@
Expand All @@ -347,7 +358,7 @@ double get_QED_2L(EvaluationContext& context)
template<class PhotonEmitter, class ExchangeField>
double AMuonVertexCorrectionFS<
PhotonEmitter, ExchangeField
>::value(const typename field_indices<EDMField>::type &indices, EvaluationContext &context)
>::value(const typename field_indices<Muon>::type &indices, const EvaluationContext &context)
{
double res = 0.0;

Expand All @@ -357,7 +368,7 @@ double AMuonVertexCorrectionFS<
PhotonEmitter
>;

constexpr auto indexBounds = FermionVertex::index_bounds;
constexpr auto indexBounds = MuonVertex::index_bounds;
for( auto indexIt = indexBounds.begin(); indexIt != indexBounds.end(); ++indexIt )
{
auto muonIndices = MuonVertex::template fieldIndices<0>(*indexIt);
Expand All @@ -367,21 +378,23 @@ double AMuonVertexCorrectionFS<

auto photonEmitterIndices = MuonVertex::template fieldIndices<2>(*indexIt);
auto exchangeFieldIndices = MuonVertex::template fieldIndices<1>(*indexIt);
auto vertex = MuonVertex::evaluate(*indexIt, context);
auto muonVertex = MuonVertex::evaluate(*indexIt, context);

auto photonEmitterMass = context.mass<PhotonEmitter>( photonEmitterIndices );
auto exchangeFieldMass = context.mass<ExchangeField>( exchangeFieldIndices );
auto muonMass = context.mass<Muon>( muonIndices );

double photonEmitterCharge = context.electric_charge<PhotonEmitter>( photonEmitterIndices );
double photonEmitterChargeCount =
- context.electric_charge<PhotonEmitter>( photonEmitterIndices ) /
context.electric_charge<Muon>( muonIndices );

const std::complex<double> zL = vertex.left();
const std::complex<double> zR = vertex.right();
const std::complex<double> zL = muonVertex.left();
const std::complex<double> zR = muonVertex.right();

const double coeffA = std::norm(zL) + std::norm(zR);
const double coeffB = (zL * std::conj(zR) + zR * std::conj(zL)).real();

const double massRatioSquared = Sqr(photonEmitterMass / exchangeMass);
const double massRatioSquared = Sqr(photonEmitterMass / exchangeFieldMass);

const double part1 = muonMass / 12.0 * coeffA * OneLoopFunctionF1C(massRatioSquared);

Expand All @@ -399,9 +412,9 @@ double AMuonVertexCorrectionFS<
}

template<class PhotonEmitter, class ExchangeField>
double AMuonVertexCorrectionFS<
double AMuonVertexCorrectionSF<
PhotonEmitter, ExchangeField
>::value(const typename field_indices<EDMField>::type &indices, EvaluationContext &context)
>::value(const typename field_indices<Muon>::type &indices, const EvaluationContext &context)
{
double res = 0.0;

Expand All @@ -411,7 +424,7 @@ double AMuonVertexCorrectionFS<
PhotonEmitter
>;

constexpr auto indexBounds = FermionVertex::index_bounds;
constexpr auto indexBounds = MuonVertex::index_bounds;
for( auto indexIt = indexBounds.begin(); indexIt != indexBounds.end(); ++indexIt )
{
auto muonIndices = MuonVertex::template fieldIndices<0>(*indexIt);
Expand All @@ -421,24 +434,26 @@ double AMuonVertexCorrectionFS<

auto photonEmitterIndices = MuonVertex::template fieldIndices<2>(*indexIt);
auto exchangeFieldIndices = MuonVertex::template fieldIndices<1>(*indexIt);
auto vertex = MuonVertex::evaluate(*indexIt, context);
auto muonVertex = MuonVertex::evaluate(*indexIt, context);

auto photonEmitterMass = context.mass<PhotonEmitter>( photonEmitterIndices );
auto exchangeFieldMass = context.mass<ExchangeField>( exchangeFieldIndices );
auto muonMass = context.mass<Muon>( muonIndices );

double photonEmitterCharge = context.electric_charge<PhotonEmitter>( photonEmitterIndices );
double photonEmitterChargeCount =
- context.electric_charge<PhotonEmitter>( photonEmitterIndices ) /
context.electric_charge<Muon>( muonIndices );

const std::complex<double> zL = muonVertex.left();
const std::complex<double> zR = muonVertex.right();

const double coeffA = std::norm(zL) + std::norm(zR);
const double coeffB = (zL * std::conj(zR) + zR * std::conj(zL)).real();

const double massRatioSquared = Sqr(exchangeMass / photonEmitterMass);
const double massRatioSquared = Sqr(exchangeFieldMass / photonEmitterMass);

const double part1 = 1.0 / 12.0 * coeffA * OneLoopFunctionF1N(massRatioSquared);
const double part2 = exchangeMass / (6.0 * muonMass) * coeffB * OneLoopFunctionF2N(massRatioSquared);
const double part2 = exchangeFieldMass / (6.0 * muonMass) * coeffB * OneLoopFunctionF2N(massRatioSquared);

const double preFactor = - oneOver16PiSqr * photonEmitterChargeCount
* muonPhysicalMass(context) * muonMass / (photonEmitterMass * photonEmitterMass);
Expand All @@ -449,4 +464,3 @@ double AMuonVertexCorrectionFS<
return res;
}
} // anonymous namespace
} // namespace flexiblesusy
8 changes: 4 additions & 4 deletions templates/a_muon2.hpp.in
View file
Expand Up @@ -19,19 +19,19 @@
// File generated at @DateAndTime@

/**
* @file @ModelName@_a_muon.hpp
* @file @ModelName@_a_muon2.hpp
*
* This file was generated at @DateAndTime@ with FlexibleSUSY
* @FlexibleSUSYVersion@ and SARAH @SARAHVersion@ .
*/

#ifndef @ModelName@_G_MUON_MINUS_2_H
#define @ModelName@_G_MUON_MINUS_2_H
#ifndef @ModelName@_A_MUON_2_H
#define @ModelName@_A_MUON_2_H

namespace flexiblesusy {
class @ModelName@_mass_eigenstates;

namespace @ModelName@_a_muon {
namespace @ModelName@_a_muon2 {
/**
* @fn calculate_a_muon
* @brief Calculates \f$a_\mu = (g-2)_\mu/2\f$ of the muon.
Expand Down

0 comments on commit 7b16e4a

Please sign in to comment.