Merge branch 'feature-2.0' into feature-2.0-semianalytic-solver

 Conflicts:
	src/standard_model_two_scale_low_scale_constraint.cpp
	src/two_scale_solver.hpp
	templates/two_scale_low_scale_constraint.hpp.in

meta/FlexibleEFTHiggsMatching.m

@@ -1,7 +1,9 @@
 BeginPackage["FlexibleEFTHiggsMatching`", {"CConversion`", "TreeMasses`", "LoopMasses`", "Constraint`", "ThresholdCorrections`", "Parameters`"}];
 
 CallSMPoleMassFunctions::usage = "";
-CalculateRunningFermionMasses::usage = "";
+CalculateRunningUpQuarkMasses::usage = "";
+CalculateRunningDownQuarkMasses::usage = "";
+CalculateRunningDownLeptonMasses::usage = "";
 FillSMFermionPoleMasses::usage = "";
 
 Begin["`Private`"];
@@ -18,7 +20,7 @@
            result
           ];
 
-CalculateRunningFermionMasses[] :=
+CalculateRunningUpQuarkMasses[] :=
     Module[{result = "", i, iStr, mq, mqFun},
            For[i = 0, i < 3, i++,
                iStr = ToString[i];
@@ -29,6 +31,11 @@
                         "sm.get_physical().MFu(" <> iStr <> ") - " <>
                         "model.get_physical().M" <> mq <> " + model.get_M" <> mqFun <> ";\n";
               ];
+           result
+          ];
+
+CalculateRunningDownQuarkMasses[] :=
+    Module[{result = "", i, iStr, mq, mqFun},
            For[i = 0, i < 3, i++,
                iStr = ToString[i];
                mq = mqFun = CConversion`RValueToCFormString[TreeMasses`GetDownQuark[i + 1, True]];
@@ -38,6 +45,11 @@
                         "sm.get_physical().MFd(" <> iStr <> ") - " <>
                         "model.get_physical().M" <> mq <> " + model.get_M" <> mqFun <> ";\n";
               ];
+           result
+          ];
+
+CalculateRunningDownLeptonMasses[] :=
+    Module[{result = "", i, iStr, mq, mqFun},
            For[i = 0, i < 3, i++,
                iStr = ToString[i];
                mq = mqFun = CConversion`RValueToCFormString[TreeMasses`GetDownLepton[i + 1, True]];

meta/FlexibleSUSY.m

@@ -1091,31 +1091,34 @@ corresponding tadpole is real or imaginary (only in models with CP
 
 WriteMatchingClass[susyScaleMatching_List, files_List] :=
     Module[{scheme = GetRenormalizationScheme[], userMatching = "",
-            gauge1Linit = "", alphaS1Lmatching = "", alphaEM1Lmatching = "",
-            setRunningFermionMasses = "", setYukawas = "",
+            alphaS1Lmatching = "", alphaEM1Lmatching = "",
+            setRunningUpQuarkMasses = "", setRunningDownQuarkMasses = "",
+            setRunningDownLeptonMasses = "", setYukawas = "",
             callAllSMPoleMassFunctions = "",
             callAllSMPoleMassFunctionsThreads = ""},
            If[FlexibleSUSY`FlexibleEFTHiggs === True,
               If[Head[susyScaleMatching] === List,
                  userMatching = Constraint`ApplyConstraints[susyScaleMatching];
                 ];
-              gauge1Linit = Parameters`CreateLocalConstRefs[
-                  ThresholdCorrections`CalculateColorCoupling[scheme] +
-                  ThresholdCorrections`CalculateElectromagneticCoupling[scheme]];
-              alphaS1Lmatching = "delta_alpha_s = alpha_s/(2.*Pi)*(" <>
-              CConversion`RValueToCFormString[ThresholdCorrections`CalculateColorCoupling[scheme]] <> ");\n";
-              alphaEM1Lmatching = "delta_alpha_em = alpha_em/(2.*Pi)*(" <>
-              CConversion`RValueToCFormString[ThresholdCorrections`CalculateElectromagneticCoupling[scheme]] <> ");\n";
-              setRunningFermionMasses           = FlexibleEFTHiggsMatching`CalculateRunningFermionMasses[];
+              alphaS1Lmatching = Parameters`CreateLocalConstRefs[ThresholdCorrections`CalculateColorCoupling[scheme]] <> "\n" <>
+                                 "delta_alpha_s += alpha_s/(2.*Pi)*(" <>
+                                 CConversion`RValueToCFormString[ThresholdCorrections`CalculateColorCoupling[scheme]] <> ");\n";
+              alphaEM1Lmatching = Parameters`CreateLocalConstRefs[ThresholdCorrections`CalculateElectromagneticCoupling[scheme]] <> "\n" <>
+                                  "delta_alpha_em += alpha_em/(2.*Pi)*(" <>
+                                  CConversion`RValueToCFormString[ThresholdCorrections`CalculateElectromagneticCoupling[scheme]] <> ");\n";
+              setRunningUpQuarkMasses           = FlexibleEFTHiggsMatching`CalculateRunningUpQuarkMasses[];
+              setRunningDownQuarkMasses         = FlexibleEFTHiggsMatching`CalculateRunningDownQuarkMasses[];
+              setRunningDownLeptonMasses        = FlexibleEFTHiggsMatching`CalculateRunningDownLeptonMasses[];
               setYukawas                        = ThresholdCorrections`SetDRbarYukawaCouplings[];
               callAllSMPoleMassFunctions        = FlexibleEFTHiggsMatching`CallSMPoleMassFunctions[FlexibleSUSY`FSEigenstates, False];
               callAllSMPoleMassFunctionsThreads = FlexibleEFTHiggsMatching`CallSMPoleMassFunctions[FlexibleSUSY`FSEigenstates, True];
              ];
            WriteOut`ReplaceInFiles[files,
-                       { "@gauge1Linit@"             -> IndentText[WrapLines[gauge1Linit]],
-                         "@alphaS1Lmatching@"        -> IndentText[WrapLines[alphaS1Lmatching]],
-                         "@alphaEM1Lmatching@"       -> IndentText[WrapLines[alphaEM1Lmatching]],
-                         "@setRunningFermionMasses@" -> IndentText[setRunningFermionMasses],
+                       { "@alphaS1Lmatching@"        -> IndentText[IndentText[WrapLines[alphaS1Lmatching]]],
+                         "@alphaEM1Lmatching@"       -> IndentText[IndentText[WrapLines[alphaEM1Lmatching]]],
+                         "@setRunningUpQuarkMasses@" -> IndentText[IndentText[setRunningUpQuarkMasses]],
+                         "@setRunningDownQuarkMasses@" -> IndentText[IndentText[setRunningDownQuarkMasses]],
+                         "@setRunningDownLeptonMasses@" -> IndentText[IndentText[setRunningDownLeptonMasses]],
                          "@setYukawas@"              -> IndentText[WrapLines[setYukawas]],
                          "@applyUserMatching@"       -> IndentText[IndentText[WrapLines[userMatching]]],
                          "@callAllSMPoleMassFunctions@" -> IndentText[callAllSMPoleMassFunctions],

src/ckm.cpp

@@ -17,6 +17,7 @@
 // ====================================================================
 
 #include "ckm.hpp"
+#include "error.hpp"
 #include "ew_input.hpp"
 #include "wrappers.hpp"
 
@@ -45,10 +46,10 @@ void CKM_parameters::reset_to_observation()
 void CKM_parameters::set_from_wolfenstein(double lambdaW, double aCkm,
                                           double rhobar, double etabar)
 {
-   assert(Abs(lambdaW) <= 1. && "Error: Wolfenstein lambda out of range!");
-   assert(Abs(aCkm)    <= 1. && "Error: Wolfenstein A parameter out of range!");
-   assert(Abs(rhobar)  <= 1. && "Error: Wolfenstein rho-bar parameter out of range!");
-   assert(Abs(etabar)  <= 1. && "Error: Wolfenstein eta-bar parameter out of range!");
+   if (Abs(lambdaW) > 1.) throw SetupError("Error: Wolfenstein lambda out of range!");
+   if (Abs(aCkm)    > 1.) throw SetupError("Error: Wolfenstein A parameter out of range!");
+   if (Abs(rhobar)  > 1.) throw SetupError("Error: Wolfenstein rho-bar parameter out of range!");
+   if (Abs(etabar)  > 1.) throw SetupError("Error: Wolfenstein eta-bar parameter out of range!");
 
    theta_12 = ArcSin(lambdaW);
    theta_23 = ArcSin(aCkm * Sqr(lambdaW));

src/eigen_utils.hpp

@@ -20,6 +20,7 @@
 #define EIGEN_UTILS_H
 
 #include <Eigen/Core>
+#include <cassert>
 #include <iomanip>
 #include <sstream>
 #include <string>

src/fixed_point_iterator.hpp

@@ -20,7 +20,6 @@
 #define FIXED_POINT_ITERATOR_H
 
 #include <iostream>
-#include <cassert>
 #include <limits>
 #include <string>
 #include <utility>
@@ -58,7 +57,7 @@ class Convergence_tester_absolute {
     * @return GSL error code (GSL_SUCCESS or GSL_CONTINUE)
     */
    int operator()(const GSL_vector& a, const GSL_vector& b) const {
-      assert(a.size() == b.size() && "Error: vectors have different size.");
+      if (a.size() != b.size()) throw SetupError("Error: vectors have different size.");
 
       const std::size_t dimension = a.size();
       double residual = 0.;
@@ -96,7 +95,7 @@ class Convergence_tester_relative {
     * @return GSL error code (GSL_SUCCESS or GSL_CONTINUE)
     */
    int operator()(const GSL_vector& a, const GSL_vector& b) const {
-      assert(a.size() == b.size() && "Error: vectors have different size.");
+      if (a.size() != b.size()) throw SetupError("Error: vectors have different size.");
 
       const std::size_t dimension = a.size();
       double rel_diff = 0.;
@@ -145,7 +144,7 @@ class Convergence_tester_tadpole {
     * @return GSL error code (GSL_SUCCESS or GSL_CONTINUE)
     */
    int operator()(const GSL_vector& a, const GSL_vector& b) const {
-      assert(a.size() == b.size() && "Error: vectors have different size.");
+      if (a.size() != b.size()) throw SetupError("Error: vectors have different size.");
 
       if (precision < 0.)
          GSL_ERROR("relative tolerance is negative", GSL_EBADTOL);

src/physical_input.cpp

@@ -17,8 +17,7 @@
 // ====================================================================
 
 #include "physical_input.hpp"
-
-#include <cassert>
+#include "error.hpp"
 
 namespace flexiblesusy {
 
@@ -63,7 +62,9 @@ void Physical_input::set(Input o, double value)
 
 void Physical_input::set(const Eigen::ArrayXd& vec)
 {
-   assert(vec.size() == values.size() && "Parameters array has wrong size");
+   if (vec.size() != values.size())
+      throw SetupError("Parameters array has wrong size");
+
    std::copy(vec.data(), vec.data() + vec.size(), values.begin());
 }
 

src/slha_io.hpp

@@ -19,7 +19,6 @@
 #ifndef SLHA_IO_H
 #define SLHA_IO_H
 
-#include <cassert>
 #include <string>
 #include <sstream>
 #include <iosfwd>
@@ -269,7 +268,7 @@ Scalar SLHA_io::convert_to(const std::string& str)
 template <class Derived>
 double SLHA_io::read_matrix(const std::string& block_name, Eigen::MatrixBase<Derived>& matrix) const
 {
-   assert(matrix.cols() > 1 && "Matrix has not more than 1 column");
+   if (matrix.cols() <= 1) throw SetupError("Matrix has less than 2 columns");
 
    SLHAea::Coll::const_iterator block =
       data.find(data.cbegin(), data.cend(), block_name);
@@ -315,7 +314,7 @@ double SLHA_io::read_matrix(const std::string& block_name, Eigen::MatrixBase<Der
 template <class Derived>
 double SLHA_io::read_vector(const std::string& block_name, Eigen::MatrixBase<Derived>& vector) const
 {
-   assert(vector.cols() == 1 && "Vector has more than 1 columns");
+   if (vector.cols() != 1) throw SetupError("Vector has more than 1 column");
 
    SLHAea::Coll::const_iterator block =
       data.find(data.cbegin(), data.cend(), block_name);

src/standard_model_two_scale_low_scale_constraint.cpp

@@ -28,21 +28,12 @@
 #include "threshold_loop_functions.hpp"
 #include "weinberg_angle.hpp"
 
-#include <cassert>
 #include <cmath>
 #include <limits>
 
 namespace flexiblesusy {
 namespace standard_model {
 
-Standard_model_low_scale_constraint<Two_scale>::Standard_model_low_scale_constraint()
-   : Single_scale_constraint()
-   , scale(0.)
-   , model(nullptr)
-   , qedqcd()
-{
-}
-
 Standard_model_low_scale_constraint<Two_scale>::Standard_model_low_scale_constraint(
    StandardModel<Two_scale>* model_, const softsusy::QedQcd& qedqcd_)
    : Single_scale_constraint()
@@ -58,8 +49,9 @@ Standard_model_low_scale_constraint<Two_scale>::~Standard_model_low_scale_constr
 
 void Standard_model_low_scale_constraint<Two_scale>::apply()
 {
-   assert(model && "Error: Standard_model_low_scale_constraint::apply():"
-          " model pointer must not be zero");
+   if (!model)
+      throw SetupError("Error: Standard_model_low_scale_constraint::apply():"
+                       " model pointer must not be zero");
 
    qedqcd.run_to(scale, 1.0e-5);
    model->calculate_DRbar_masses();

src/standard_model_two_scale_low_scale_constraint.hpp

@@ -37,7 +37,7 @@ class StandardModel;
 template<>
 class Standard_model_low_scale_constraint<Two_scale> : public Single_scale_constraint {
 public:
-   Standard_model_low_scale_constraint();
+   Standard_model_low_scale_constraint() = default;
    Standard_model_low_scale_constraint(StandardModel<Two_scale>*, const softsusy::QedQcd&);
    virtual ~Standard_model_low_scale_constraint();
    virtual void apply() override;
@@ -52,9 +52,9 @@ class Standard_model_low_scale_constraint<Two_scale> : public Single_scale_const
    void set_threshold_corrections_loop_order(unsigned); ///< threshold corrections loop order
 
 private:
-   double scale;
-   StandardModel<Two_scale>* model;
-   softsusy::QedQcd qedqcd;
+   double scale{0.};
+   StandardModel<Two_scale>* model{nullptr};
+   softsusy::QedQcd qedqcd{};
 };
 
 } // namespace standard_model

src/two_scale_solver.hpp

@@ -61,8 +61,6 @@ class Two_scale_running_precision;
 template<>
 class RGFlow<Two_scale> {
 public:
-   /// Create empty two scale solver.
-   /// The RG running precision is set to the default value 0.001.
    RGFlow() = default;
    RGFlow(const RGFlow&) = delete;
    RGFlow(RGFlow&&) = delete;
@@ -145,12 +143,11 @@ class RGFlow<Two_scale> {
    void clear_problems();              ///< clear model problems
    unsigned int get_max_iterations() const; ///< returns max. number of iterations
    Model* get_model(double) const;     ///< returns model at given scale
-   void initial_guess();               ///< initial guess
    double get_precision();             ///< returns running precision
-   void update_running_precision();    ///< update the RG running precision
-   std::vector<std::shared_ptr<Slider> > sort_sliders() const; ///< sort the sliders w.r.t. to scale
-
+   void initial_guess();               ///< initial guess
    void run_sliders();                 ///< run all sliders
+   std::vector<std::shared_ptr<Slider> > sort_sliders() const; ///< sort the sliders w.r.t. to scale
+   void update_running_precision();    ///< update the RG running precision
 };
 
 }

src/wrappers.hpp

@@ -34,6 +34,7 @@
 
 #include "dilog.hpp"
 #include "eigen_tensor.hpp"
+#include "error.hpp"
 #include "if.hpp"
 #include "sum.hpp"
 #include "which.hpp"
@@ -348,7 +349,8 @@ double MaxRelDiff(const Eigen::MatrixBase<Derived>& a,
 {
    typename Eigen::MatrixBase<Derived>::PlainObject sumTol(a.rows());
 
-   assert(a.rows() == b.rows());
+   if (a.rows() != b.rows())
+      throw SetupError("MaxRelDiff: vectors have different size!");
 
    for (int i = 0; i < a.rows(); i++)
       sumTol(i) = MaxRelDiff(a(i), b(i));
@@ -399,7 +401,7 @@ template <typename T>
 T PolyLog(int n, T z) {
    if (n == 2)
       return gm2calc::dilog(z);
-   assert(false && "PolyLog(n!=2) not implemented");
+   throw SetupError("PolyLog(n!=2) not implemented");
 }
 
 template <typename Base, typename Exponent>