remove QedQcd member from Standard_model class

as it is used only for initialization.

This change avoids bugs like the ones fixed in commit 20f169f by
design.  In addition the Standard_model class is now a bit smaller.

In addition the calculation of the Standard_model effective couplings
was optimized as in commit b533d67 .

src/standard_model.cpp

@@ -26,6 +26,7 @@
 #include "eigen_utils.hpp"
 #include "wrappers.hpp"
 #include "linalg2.hpp"
+#include "lowe.h"
 #include "numerics2.hpp"
 #include "logger.hpp"
 #include "error.hpp"
@@ -163,7 +164,6 @@ Standard_model::Standard_model()
    , physical()
    , problems(standard_model_info::particle_names)
    , two_loop_corrections()
-   , qedqcd()
    , input()
    , g1(0), g2(0), g3(0), Lambdax(0), Yu(Eigen::Matrix<double,3,3>::Zero()), Yd
    (Eigen::Matrix<double,3,3>::Zero()), Ye(Eigen::Matrix<double,3,3>::Zero())
@@ -198,7 +198,6 @@ Standard_model::Standard_model(double scale_, double loops_, double thresholds_
    , physical()
    , problems(standard_model_info::particle_names)
    , two_loop_corrections()
-   , qedqcd()
    , input()
    , g1(g1_), g2(g2_), g3(g3_), Lambdax(Lambdax_), Yu(Yu_), Yd(Yd_), Ye(Ye_)
    , mu2(mu2_), v(v_)
@@ -869,13 +868,14 @@ std::string Standard_model::name() const
    return "Standard model";
 }
 
-void Standard_model::initialise_from_input()
+void Standard_model::initialise_from_input(const softsusy::QedQcd& qedqcd_)
 {
    const double scale = get_scale();
+   auto qedqcd = qedqcd_;
 
    // initial guess
    qedqcd.to(qedqcd.displayPoleMZ());
-   initial_guess_for_parameters();
+   initial_guess_for_parameters(qedqcd);
    run_to(qedqcd.displayPoleMZ());
 
    // determine Standard model parameters iteratively
@@ -907,7 +907,7 @@ void Standard_model::initialise_from_input()
       const double alpha_em_drbar = alpha_em / (1.0 - delta_alpha_em);
       const double alpha_s_drbar  = alpha_s / (1.0 - delta_alpha_s);
       const double e_drbar        = Sqrt(4.0 * Pi * alpha_em_drbar);
-      double theta_w_drbar        = calculate_theta_w(alpha_em_drbar);
+      double theta_w_drbar        = calculate_theta_w(qedqcd, alpha_em_drbar);
 
       if (IsFinite(theta_w_drbar)) {
          problems.unflag_non_perturbative_parameter(
@@ -920,9 +920,9 @@ void Standard_model::initialise_from_input()
 
       v = Re((2 * mz_run) / Sqrt(0.6 * Sqr(g1) + Sqr(g2)));
 
-      calculate_Yu_DRbar();
-      calculate_Yd_DRbar();
-      calculate_Ye_DRbar();
+      calculate_Yu_DRbar(qedqcd);
+      calculate_Yd_DRbar(qedqcd);
+      calculate_Ye_DRbar(qedqcd);
       calculate_Lambdax_DRbar();
 
       solve_ewsb();
@@ -931,7 +931,8 @@ void Standard_model::initialise_from_input()
       g2 = e_drbar * Csc(theta_w_drbar);
       g3 = 3.5449077018110318 * Sqrt(alpha_s_drbar);
 
-      recalculate_mw_pole();
+      if (get_thresholds())
+         qedqcd.setPoleMW(recalculate_mw_pole(qedqcd.displayPoleMW()));
 
       converged = check_convergence(old);
       old = *this;
@@ -942,7 +943,7 @@ void Standard_model::initialise_from_input()
       run_to(scale);
 }
 
-void Standard_model::initial_guess_for_parameters()
+void Standard_model::initial_guess_for_parameters(const softsusy::QedQcd& qedqcd)
 {
    const double MZ = qedqcd.displayPoleMZ();
    const double MW = qedqcd.displayPoleMW();
@@ -1017,7 +1018,7 @@ double Standard_model::calculate_delta_alpha_s(double alphaS) const
 
 }
 
-double Standard_model::calculate_theta_w(double alpha_em_drbar)
+double Standard_model::calculate_theta_w(const softsusy::QedQcd& qedqcd, double alpha_em_drbar)
 {
    double theta_w = 0.;
 
@@ -1093,7 +1094,7 @@ double Standard_model::calculate_theta_w(double alpha_em_drbar)
    return theta_w;
 }
 
-void Standard_model::calculate_Yu_DRbar()
+void Standard_model::calculate_Yu_DRbar(const softsusy::QedQcd& qedqcd)
 {
    Eigen::Matrix<double,3,3> upQuarksDRbar(Eigen::Matrix<double,3,3>::Zero());
    upQuarksDRbar(0,0)      = qedqcd.displayMass(softsusy::mUp);
@@ -1107,7 +1108,7 @@ void Standard_model::calculate_Yu_DRbar()
    Yu = -((1.4142135623730951*upQuarksDRbar)/v).transpose();
 }
 
-void Standard_model::calculate_Yd_DRbar()
+void Standard_model::calculate_Yd_DRbar(const softsusy::QedQcd& qedqcd)
 {
    Eigen::Matrix<double,3,3> downQuarksDRbar(Eigen::Matrix<double,3,3>::Zero());
    downQuarksDRbar(0,0)   = qedqcd.displayMass(softsusy::mDown);
@@ -1121,7 +1122,7 @@ void Standard_model::calculate_Yd_DRbar()
    Yd = ((1.4142135623730951*downQuarksDRbar)/v).transpose();
 }
 
-void Standard_model::calculate_Ye_DRbar()
+void Standard_model::calculate_Ye_DRbar(const softsusy::QedQcd& qedqcd)
 {
    Eigen::Matrix<double,3,3> downLeptonsDRbar(Eigen::Matrix<double,3,3>::Zero());
    downLeptonsDRbar(0,0) = qedqcd.displayPoleMel();
@@ -1148,12 +1149,6 @@ void Standard_model::calculate_Lambdax_DRbar()
    Lambdax = Sqr(higgsDRbar) / Sqr(v);
 }
 
-void Standard_model::recalculate_mw_pole()
-{
-   if (get_thresholds())
-      qedqcd.setPoleMW(recalculate_mw_pole(qedqcd.displayPoleMW()));
-}
-
 double Standard_model::recalculate_mw_pole(double mw_pole)
 {
    calculate_MVWp();

src/standard_model.hpp

@@ -32,7 +32,6 @@
 #include "error.hpp"
 #include "problems.hpp"
 #include "config.h"
-#include "lowe.h"
 #include "physical_input.hpp"
 
 #include <iosfwd>
@@ -41,6 +40,10 @@
 #include <gsl/gsl_vector.h>
 #include <Eigen/Core>
 
+namespace softsusy {
+   class QedQcd;
+}
+
 namespace flexiblesusy {
 
 class EWSB_solver;
@@ -151,14 +154,11 @@ class Standard_model : public Beta_function {
    void set_precision(double);
    double get_precision() const;
 
-   void set_low_energy_data(const softsusy::QedQcd& qedqcd_) { qedqcd = qedqcd_; }
-   const softsusy::QedQcd& get_low_energy_data() const { return qedqcd; }
-   softsusy::QedQcd& get_low_energy_data() { return qedqcd; }
    void set_physical_input(const Physical_input& input_) { input = input_; }
    const Physical_input& get_physical_input() const { return input; }
    Physical_input& get_physical_input() { return input; }
 
-   void initialise_from_input();
+   void initialise_from_input(const softsusy::QedQcd&);
 
    void set_g1(double g1_) { g1 = g1_; }
    void set_g2(double g2_) { g2 = g2_; }
@@ -485,11 +485,10 @@ class Standard_model : public Beta_function {
    double calculate_delta_alpha_em(double alphaEm) const;
    double calculate_delta_alpha_s(double alphaS) const;
    void calculate_Lambdax_DRbar();
-   double calculate_theta_w(double alpha_em_drbar);
-   void calculate_Yu_DRbar();
-   void calculate_Yd_DRbar();
-   void calculate_Ye_DRbar();
-   void recalculate_mw_pole();
+   double calculate_theta_w(const softsusy::QedQcd&, double alpha_em_drbar);
+   void calculate_Yu_DRbar(const softsusy::QedQcd&);
+   void calculate_Yd_DRbar(const softsusy::QedQcd&);
+   void calculate_Ye_DRbar(const softsusy::QedQcd&);
    double recalculate_mw_pole(double);
 
 protected:
@@ -575,7 +574,6 @@ class Standard_model : public Beta_function {
    Standard_model_physical physical; ///< contains the pole masses and mixings
    Problems<standard_model_info::NUMBER_OF_PARTICLES> problems;
    Two_loop_corrections two_loop_corrections; ///< used 2-loop corrections
-   softsusy::QedQcd qedqcd;
    Physical_input input;
 
    int solve_ewsb_iteratively();
@@ -588,7 +586,7 @@ class Standard_model : public Beta_function {
    static int tadpole_equations(const gsl_vector*, void*, gsl_vector*);
    void copy_DRbar_masses_to_pole_masses();
 
-   void initial_guess_for_parameters();
+   void initial_guess_for_parameters(const softsusy::QedQcd&);
    bool check_convergence(const Standard_model& old) const;
 
    // Passarino-Veltman loop functions

src/standard_model_effective_couplings.cpp

@@ -53,16 +53,17 @@ Standard_model_effective_couplings::~Standard_model_effective_couplings()
 
 void Standard_model_effective_couplings::calculate_effective_couplings()
 {
+   const standard_model::Standard_model sm(initialise_SM());
    const double scale = model.get_scale();
    const Eigen::ArrayXd saved_parameters(model.get());
 
    const double saved_mt = PHYSICAL(MFu(2));
    PHYSICAL(MFu(2)) = qedqcd.displayPoleMt();
 
    const auto Mhh = PHYSICAL(Mhh);
-   run_SM_strong_coupling_to(0.5 * Mhh);
+   run_SM_strong_coupling_to(sm, 0.5 * Mhh);
    calculate_eff_CphhVPVP();
-   run_SM_strong_coupling_to(Mhh);
+   run_SM_strong_coupling_to(sm, Mhh);
    calculate_eff_CphhVGVG();
 
    PHYSICAL(MFu(2)) = saved_mt;
@@ -89,22 +90,23 @@ void Standard_model_effective_couplings::copy_mixing_matrices_from_model()
 
 }
 
-void Standard_model_effective_couplings::run_SM_strong_coupling_to(double m)
+standard_model::Standard_model Standard_model_effective_couplings::initialise_SM() const
 {
-   using namespace standard_model;
-
-   Standard_model sm;
+   standard_model::Standard_model sm;
 
    sm.set_loops(2);
    sm.set_thresholds(2);
-   sm.set_low_energy_data(qedqcd);
    sm.set_physical_input(physical_input);
 
-   sm.initialise_from_input();
-   sm.run_to(m);
+   sm.initialise_from_input(qedqcd);
 
-   model.set_g3(sm.get_g3());
+   return sm;
+}
 
+void Standard_model_effective_couplings::run_SM_strong_coupling_to(standard_model::Standard_model sm, double m)
+{
+   sm.run_to(m);
+   model.set_g3(sm.get_g3());
 }
 
 std::complex<double> Standard_model_effective_couplings::scalar_scalar_qcd_factor(double m_decay, double m_loop) const

src/standard_model_effective_couplings.hpp

@@ -78,7 +78,8 @@ class Standard_model_effective_couplings {
 
    void copy_mixing_matrices_from_model();
 
-   void run_SM_strong_coupling_to(double m);
+   standard_model::Standard_model initialise_SM() const;
+   void run_SM_strong_coupling_to(standard_model::Standard_model, double m);
 
    // higher order corrections to the amplitudes for
    // effective coupling to photons

src/standard_model_two_scale_low_scale_constraint.cpp

@@ -62,7 +62,6 @@ void Standard_model_low_scale_constraint<Two_scale>::apply()
           " model pointer must not be zero");
 
    qedqcd.runto(scale, 1.0e-5);
-   model->set_low_energy_data(qedqcd);
    model->calculate_DRbar_masses();
 
    const double alpha_em = qedqcd.displayAlpha(softsusy::ALPHA);
@@ -84,7 +83,7 @@ void Standard_model_low_scale_constraint<Two_scale>::apply()
    const double g2 = model->get_g2();
    const double mZ = model->get_thresholds() ?
       model->calculate_MVZ_DRbar(mz_pole) : mz_pole;
-   double theta_w = model->calculate_theta_w(alpha_em_drbar);
+   double theta_w = model->calculate_theta_w(qedqcd, alpha_em_drbar);
 
    if (IsFinite(theta_w)) {
       model->get_problems().unflag_non_perturbative_parameter(
@@ -96,9 +95,9 @@ void Standard_model_low_scale_constraint<Two_scale>::apply()
    }
 
    model->set_v(Re((2*mZ)/Sqrt(0.6*Sqr(g1) + Sqr(g2))));
-   model->calculate_Yu_DRbar();
-   model->calculate_Yd_DRbar();
-   model->calculate_Ye_DRbar();
+   model->calculate_Yu_DRbar(qedqcd);
+   model->calculate_Yd_DRbar(qedqcd);
+   model->calculate_Ye_DRbar(qedqcd);
    model->set_g1(1.2909944487358056*e_drbar*Sec(theta_w));
    model->set_g2(e_drbar*Csc(theta_w));
    model->set_g3(3.5449077018110318*Sqrt(alpha_s_drbar));

templates/effective_couplings.cpp.in

@@ -68,10 +68,9 @@ standard_model::Standard_model @ModelName@_effective_couplings::initialise_SM()
 
    sm.set_loops(2);
    sm.set_thresholds(2);
-   sm.set_low_energy_data(qedqcd);
    sm.set_physical_input(physical_input);
 
-   sm.initialise_from_input();
+   sm.initialise_from_input(qedqcd);
 
    return sm;
 }

test/test_MSSMtower.sh

@@ -210,7 +210,7 @@ EOF
     done
 }
 
-scan "$start" "$stop" "$steps" | tee "$scan_data"
+# scan "$start" "$stop" "$steps" | tee "$scan_data"
 
 cat <<EOF | gnuplot
 set terminal pdf enhanced