Skip to content

Commit

Permalink
Merge pull request #20 from FlexibleSUSY/feature-2.0-amu-fix
Browse files Browse the repository at this point in the history
Feature 2.0 amu fix
  • Loading branch information
Expander committed Sep 11, 2017
2 parents 82dc9f6 + f229034 commit a523025
Show file tree
Hide file tree
Showing 6 changed files with 170 additions and 26 deletions.
37 changes: 19 additions & 18 deletions meta/GMuonMinus2.m
Expand Up @@ -357,20 +357,12 @@ If you add new kinds of vertices (e.g for new diagram types):
(* Use of memoization gives ~30% speedup for the MSSM! *)
memoizedVertices = {};
MemoizingVertex[particles_List, options : OptionsPattern[SARAH`Vertex]] :=
Module[{memo, ordering, orderedParticles},
(* First we sort the particles *)
ordering = Ordering[particles];
orderedParticles = particles[[ordering]];

memo = Select[memoizedVertices, MatchesMemoizedVertex[orderedParticles], 1];
If[memo =!= {}, memo = memo[[1]],
Module[{memo = Select[memoizedVertices, MatchesMemoizedVertex[particles], 1]},
If[memo =!= {}, memo[[1]],
(* Create a new entry *)
memo = SARAH`Vertex[orderedParticles, options];
AppendTo[memoizedVertices, memo];];

(* Now return the particles to their original order *)
memo = OrderVertex[memo, Ordering[ordering]];
Return[memo]];
memo = SARAH`Vertex[particles, options];
AppendTo[memoizedVertices, memo];
memo]];

MatchesMemoizedVertex[particles_List][vertex_] := MatchQ[particles, Vertices`StripFieldIndices /@ vertex[[1]]];

Expand Down Expand Up @@ -567,7 +559,7 @@ If you add new kinds of vertices (e.g for new diagram types):
">\n" <>
"{\n" <>
IndentText @
("static const bool is_permutation = true;\n" <>
("static const int permutation = " <> ToString[PermutationFactor[particles, ordering]] <> ";\n" <>
"typedef " <> orderedVertexFunction <> " orig_type;\n" <>
"typedef boost::mpl::vector_c<unsigned, " <>
StringJoin @ Riffle[ToString /@ (Ordering[ordering] - 1), ", "] <>
Expand All @@ -579,6 +571,15 @@ If you add new kinds of vertices (e.g for new diagram types):
Return[{prototypes, definitions}];
];

PermutationFactor[particles_List, ordering_List] := Module[{
particleTypes = SARAH`getType[#, False, FlexibleSUSY`FSEigenstates]& /@
particles
},
Which[SARAH`VType @@ particleTypes =!= FFV, 1,
OrderedQ @ Select[ordering, particleTypes[[#]] === F &], 1,
True, -1]
];

(* Creates local declarations of field indices, whose values are taken
from the elements of `arrayName'.
*)
Expand Down Expand Up @@ -628,7 +629,7 @@ If you add new kinds of vertices (e.g for new diagram types):
Parameters`CreateLocalConstRefs[expr] <> "\n" <>
"const " <> GetComplexScalarCType[] <> " result = " <>
Parameters`ExpressionToString[expr] <> ";\n\n" <>
"return vertex_type(result);",
"return vertex_type(result * permutationFactor);",

"LeftAndRightComponentedVertex",
exprL = Vertices`SortCp @ SARAH`Cp[Sequence @@ indexedParticles][SARAH`PL] /. vertexRules;
Expand All @@ -641,7 +642,7 @@ If you add new kinds of vertices (e.g for new diagram types):
Parameters`ExpressionToString[exprL] <> ";\n\n" <>
"const " <> GetComplexScalarCType[] <> " right = " <>
Parameters`ExpressionToString[exprR] <> ";\n\n" <>
"return vertex_type(left, right);"];
"return vertex_type(left * permutationFactor, right * permutationFactor);"];

sarahParticles = SARAH`getParticleName /@ particles;
particleInfo = Flatten[(Cases[SARAH`Particles[FlexibleSUSY`FSEigenstates], {#, ___}] &) /@
Expand Down Expand Up @@ -699,7 +700,7 @@ If you add new kinds of vertices (e.g for new diagram types):
prototype = ("template<> struct " <> dataClassName <> "\n" <>
"{\n" <>
IndentText @
("static const bool is_permutation = false;\n" <>
("static const int permutation = 0;\n" <>
"typedef IndexBounds<" <> ToString @ NumberOfIndices[parsedVertex] <> "> index_bounds;\n" <>
"typedef " <> VertexClassName[parsedVertex] <> " vertex_type;\n" <>
"typedef boost::mpl::vector_c<unsigned, " <>
Expand All @@ -723,7 +724,7 @@ If you add new kinds of vertices (e.g for new diagram types):
);
];
definition = ("template<> template<> " <> functionClassName <> "::vertex_type\n" <>
functionClassName <> "::vertex(const indices_type &indices, EvaluationContext &context)\n" <>
functionClassName <> "::vertex(const indices_type &indices, EvaluationContext &context, int permutationFactor)\n" <>
"{\n" <>
IndentText @ VertexFunctionBody[parsedVertex] <> "\n" <>
"}");
Expand Down
5 changes: 5 additions & 0 deletions model_files/MRSSM2/FlexibleSUSY.m.in
Expand Up @@ -87,6 +87,11 @@ InitialGuessAtLowScale = {
{Ye, Automatic}
};

ExtraSLHAOutputBlocks = {
{FlexibleSUSYLowEnergy,
{{21, FlexibleSUSYObservable`aMuon} } }
};

DefaultPoleMassPrecision = HighPrecision;
HighPoleMassPrecision = {hh, Ah, Hpm};
MediumPoleMassPrecision = {};
Expand Down
20 changes: 12 additions & 8 deletions templates/a_muon.cpp.in
Expand Up @@ -296,15 +296,15 @@ public:
template<class ...P> struct VertexFunctionData;

/**
* @class VertexFunctionHelper<bool, P...>
* @class VertexFunctionHelper<int, P...>
* @brief A helper template that is used by VertexFunction<P...>
*
* You should never have to use this class directly. It is just a
* helper class.
*/
template<bool, class ...P> class VertexFunctionHelper;
template<int, class ...P> class VertexFunctionHelper;

template<class ...P> class VertexFunctionHelper<false, P...>
template<class ...P> class VertexFunctionHelper<0, P...>
: public VertexFunctionData<P...> {
private:
typedef VertexFunctionData<P...> Base;
Expand All @@ -326,9 +326,11 @@ public:
auto end = indices.begin() + boost::mpl::at_c<particleIndexStart, particleIndex+1>::type::value;
return std::vector<unsigned>(begin, end);
}

static constexpr int permutationFactor = 1;
};

template<class ...P> class VertexFunctionHelper<true, P...>
template<int permutation, class ...P> class VertexFunctionHelper
: public VertexFunctionData<P...>::orig_type {
private:
typedef typename VertexFunctionData<P...>::orig_type orig_type;
Expand All @@ -345,6 +347,8 @@ public:
boost::mpl::at_c<particlePermutation, particleIndex>::type::value
> (indices);
}

static constexpr int permutationFactor = permutation;
};

/**
Expand All @@ -359,9 +363,9 @@ public:
* possible particle permutations.
*/
template<class ...P> class VertexFunction
: public VertexFunctionHelper<VertexFunctionData<P...>::is_permutation, P...> {
: public VertexFunctionHelper<VertexFunctionData<P...>::permutation, P...> {
private:
typedef VertexFunctionHelper<VertexFunctionData<P...>::is_permutation, P...> Base;
typedef VertexFunctionHelper<VertexFunctionData<P...>::permutation, P...> Base;
public:
typedef typename Base::indices_type indices_type;
typedef typename Base::index_bounds index_bounds;
Expand Down Expand Up @@ -395,9 +399,9 @@ public:
* For canonical P... the definitions have to be supplied separately.
*/
template<class Q = Base>
static vertex_type vertex(const indices_type &indices, EvaluationContext &context)
static vertex_type vertex(const indices_type &indices, EvaluationContext &context, int permutationFactor = 1)
{
return Q::vertex(indices, context);
return Q::vertex(indices, context, Base::permutationFactor);
}
};

Expand Down
7 changes: 7 additions & 0 deletions test/module.mk
Expand Up @@ -180,6 +180,11 @@ TEST_SRC += \
$(DIR)/test_CMSSM_database.cpp
endif

ifeq ($(WITH_MRSSM2),yes)
TEST_SRC += \
$(DIR)/test_MRSSM2_gmm2.cpp
endif

endif # ifneq ($(findstring two_scale,$(SOLVERS)),)

ifneq ($(findstring semi_analytic,$(SOLVERS)),)
Expand Down Expand Up @@ -800,6 +805,8 @@ $(DIR)/test_sfermions.x: $(LIBSoftsusyMSSM) $(LIBCMSSM) $(LIBFLEXI) $(LIBTEST) $
$(DIR)/test_CMSSM_database.x: $(DIR)/test_CMSSM_database.o $(LIBCMSSM) $(LIBFLEXI) $(LIBTEST) $(filter-out -%,$(LOOPFUNCLIBS))
$(CXX) $(CXXFLAGS) $(CPPFLAGS) -o $@ $(call abspathx,$^) $(BOOSTTESTLIBS) $(BOOSTTHREADLIBS) $(GSLLIBS) $(FLIBS) $(SQLITELIBS) $(THREADLIBS)

$(DIR)/test_MRSSM2_gmm2.x: $(LIBMRSSM2) $(LIBFLEXI) $(filter-out -%,$(LOOPFUNCLIBS))

$(DIR)/test_CMSSM_model.x: $(LIBSoftsusyMSSM) $(LIBCMSSM) $(LIBFLEXI) $(LIBTEST) $(filter-out -%,$(LOOPFUNCLIBS))

$(DIR)/test_CMSSM_info.x: $(LIBCMSSM) $(LIBFLEXI) $(filter-out -%,$(LOOPFUNCLIBS))
Expand Down
45 changes: 45 additions & 0 deletions test/test_MRSSM2.hpp
@@ -0,0 +1,45 @@
// ====================================================================
// 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
// <http://www.gnu.org/licenses/>.
// ====================================================================

#ifndef TEST_MRSSM2_H
#define TEST_MRSSM2_H

#include <tuple>

#include "lowe.h"
#include "MRSSM2_two_scale_spectrum_generator.hpp"

using namespace flexiblesusy;

MRSSM2_slha<MRSSM2<Two_scale>>
setup_MRSSM2(const MRSSM2_input_parameters& input)
{
softsusy::QedQcd qedqcd;

Spectrum_generator_settings settings;
settings.set(Spectrum_generator_settings::calculate_sm_masses, 0);
settings.set(Spectrum_generator_settings::calculate_bsm_masses, 0);

MRSSM2_spectrum_generator<Two_scale> spectrum_generator;
spectrum_generator.set_settings(settings);
spectrum_generator.run(qedqcd, input);

return std::get<0>(spectrum_generator.get_models_slha());
}

#endif
82 changes: 82 additions & 0 deletions test/test_MRSSM2_gmm2.cpp
@@ -0,0 +1,82 @@
// ====================================================================
// 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
// <http://www.gnu.org/licenses/>.
// ====================================================================

#define BOOST_TEST_DYN_LINK
#define BOOST_TEST_MODULE test_MRSSM2_gmm2

#include <boost/test/unit_test.hpp>

#include "test_MRSSM2.hpp"

#include "wrappers.hpp"
#include "MRSSM2_a_muon.hpp"

using namespace flexiblesusy;

BOOST_AUTO_TEST_CASE( test_amu )
{
MRSSM2_input_parameters input;

// chargino dominance
input.TanBeta = 10;
input.Ms = 1000;
input.LamTDInput = -1.0;
input.LamTUInput = -1.0;
input.LamSDInput = 1.1;
input.LamSUInput = -1.1;
input.MuDInput = 400;
input.MuUInput = 400;
input.BMuInput = Sqr(300);
input.mq2Input << Sqr(1000), 0, 0,
0, Sqr(1000), 0,
0, 0, Sqr(1000);
input.ml2Input << Sqr(500), 0, 0,
0, Sqr(500), 0,
0, 0, Sqr(500);
input.md2Input << Sqr(1000), 0, 0,
0, Sqr(1000), 0,
0, 0, Sqr(1000);
input.mu2Input << Sqr(1000), 0, 0,
0, Sqr(1000), 0,
0, 0, Sqr(1000);
input.me2Input << Sqr(500), 0, 0,
0, Sqr(500), 0,
0, 0, Sqr(500);
input.mS2Input = Sqr(2000);
input.mT2Input = Sqr(3000);
input.moc2Input = Sqr(1000);
input.mRd2Input = Sqr(700);
input.mRu2Input = Sqr(1000);
input.MDBSInput = 1000;
input.MDWBTInput = 500;
input.MDGocInput = 1500;

MRSSM2_slha<MRSSM2<Two_scale>> m = setup_MRSSM2(input);

auto amu = MRSSM2_a_muon::calculate_a_muon(m);
BOOST_CHECK_CLOSE(amu, -8.30e-11, 1.0);

// neutralino dominance
input.ml2Input << Sqr(8000), 0, 0,
0, Sqr(8000), 0,
0, 0, Sqr(8000);
m = setup_MRSSM2(input);

amu = MRSSM2_a_muon::calculate_a_muon(m);
BOOST_CHECK_CLOSE(amu, 6.33e-12, 1.0);
}

0 comments on commit a523025

Please sign in to comment.