From 3a585aaf587c3a07dec675c6f8466a3033167ee9 Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Wed, 20 Apr 2016 17:42:20 +0200 Subject: [PATCH] make to() re-callable by adding fixed input parameters --- src/lowe.cpp | 79 ++++++++++++++++++++++++++++++++++------------ src/lowe.h | 16 ++++++++++ src/slha_io.cpp | 5 +++ test/test_lowe.cpp | 21 ++++++++++++ 4 files changed, 100 insertions(+), 21 deletions(-) diff --git a/src/lowe.cpp b/src/lowe.cpp index cc78136ea..eb808892f 100644 --- a/src/lowe.cpp +++ b/src/lowe.cpp @@ -9,7 +9,6 @@ #include "lowe.h" #include "ew_input.hpp" #include "wrappers.hpp" -#include "error.hpp" namespace softsusy { @@ -17,7 +16,7 @@ const char* QedQcd_input_parmeter_names[NUMBER_OF_LOW_ENERGY_INPUT_PARAMETERS] = "alpha_em_MSbar_at_MZ", "GFermi", "alpha_s_MSbar_at_MZ", "MZ_pole", "mb_mb", "MT_pole", "MTau_pole", "MMuon_pole", "MElectron_pole", "Mv3_pole", "MW_pole", "ME_pole", "Mv1_pole", "MM_pole", "Mv2_pole", "MD_2GeV", "MU_2GeV", "MS_2GeV", - "MC_2GeV" }; + "MC_2GeV", "mc_mc", "MU_2GeV_Input", "MD_2GeV_Input", "MS_2GeV_Input" }; /// external object temp used to get objects into external routines, however: /// don't use it! @@ -33,6 +32,10 @@ QedQcd::QedQcd() , gfermi(flexiblesusy::Electroweak_constants::gfermi) , ckm() , pmns() + , mcMc(MCHARM) + , mu2GeV(MUP) + , md2GeV(MDOWN) + , ms2GeV(MSTRANGE) { setPars(11); // Default object: 1998 PDB defined in 'def.h' @@ -62,6 +65,10 @@ const QedQcd & QedQcd::operator=(const QedQcd & m) { mnu = m.mnu; ckm = m.ckm; pmns = m.pmns; + mcMc = m.mcMc; + mu2GeV = m.mu2GeV; + md2GeV = m.md2GeV; + ms2GeV = m.ms2GeV; setLoops(m.displayLoops()); setThresholds(m.displayThresholds()); setMu(m.displayMu()); @@ -422,42 +429,58 @@ void QedQcd::to(double scale, double precision_goal, unsigned max_iterations) { bool converged = false; Eigen::ArrayXd qedqcd_old(display_input()), qedqcd_new(display_input()); + runto(displayPoleMZ(), precision_goal); + // save user-defined values const QedQcd saved(*this); while (!converged && it < max_iterations) { // set alpha_i(MZ) error = runto(saved.displayPoleMZ(), precision_goal); - if (error) - throw flexiblesusy::NonPerturbativeRunningError(saved.displayPoleMZ()); + if (error) { + throw std::string("Error: Non-perturbative running to MZ = ") + + flexiblesusy::ToString(saved.displayPoleMZ()) + + " during determination of the SM(5) parameters."; + } setAlpha(ALPHA, saved.displayAlpha(ALPHA)); setAlpha(ALPHAS, saved.displayAlpha(ALPHAS)); // set mb(mb) error = runto(saved.displayMbMb(), precision_goal); - if (error) - throw flexiblesusy::NonPerturbativeRunningError(saved.displayMbMb()); + if (error) { + throw std::string("Error: Non-perturbative running to mb(mb) = ") + + flexiblesusy::ToString(saved.displayMbMb()) + + " during determination of the SM(5) parameters."; + } setMass(mBottom, saved.displayMbMb()); setPoleMb(extractPoleMb(displayAlpha(ALPHAS))); // set mc(mc) error = runto(saved.displayMass(mCharm), precision_goal); - if (error) - throw flexiblesusy::NonPerturbativeRunningError(saved.displayMass(mCharm)); - setMass(mCharm, saved.displayMass(mCharm)); + if (error) { + throw std::string("Error: Non-perturbative running to mc(mc) = ") + + flexiblesusy::ToString(saved.displayMass(mCharm)) + + " during determination of the SM(5) parameters."; + } + setMass(mCharm, saved.displayMcMc()); // set mu, md, ms at 2 GeV error = runto(2.0, precision_goal); - if (error) - throw flexiblesusy::NonPerturbativeRunningError(2.0); - setMass(mUp, saved.displayMass(mUp)); - setMass(mDown, saved.displayMass(mDown)); - setMass(mStrange, saved.displayMass(mStrange)); + if (error) { + throw std::string("Error: Non-perturbative running to 2 GeV" + " during determination of the SM(5) parameters."); + } + setMass(mUp, saved.displayMu2GeV()); + setMass(mDown, saved.displayMd2GeV()); + setMass(mStrange, saved.displayMs2GeV()); // check convergence error = runto(scale, precision_goal); - if (error) - throw flexiblesusy::NonPerturbativeRunningError(scale); + if (error) { + throw std::string("Error: Non-perturbative running to Q = ") + + flexiblesusy::ToString(scale) + + " during determination of the SM(5) parameters."; + } qedqcd_new = display_input(); converged = flexiblesusy::MaxRelDiff(qedqcd_old, qedqcd_new) < precision_goal; @@ -469,14 +492,20 @@ void QedQcd::to(double scale, double precision_goal, unsigned max_iterations) { // set alpha_i(MZ) on last time error = runto(saved.displayPoleMZ(), precision_goal); - if (error) - throw flexiblesusy::NonPerturbativeRunningError(saved.displayPoleMZ()); + if (error) { + throw std::string("Error: Non-perturbative running to MZ = ") + + flexiblesusy::ToString(saved.displayPoleMZ()) + + " during determination of the SM(5) parameters."; + } setAlpha(ALPHA, saved.displayAlpha(ALPHA)); setAlpha(ALPHAS, saved.displayAlpha(ALPHAS)); error = runto(scale, precision_goal); - if (error) - throw flexiblesusy::NonPerturbativeRunningError(scale); + if (error) { + throw std::string("Error: Non-perturbative running to Q = ") + + flexiblesusy::ToString(scale) + + " during determination of the SM(5) parameters."; + } if (!converged && max_iterations > 0) { std::ostringstream ostr; @@ -706,11 +735,15 @@ void QedQcd::set_input(const Eigen::ArrayXd& pars) mf(1) = pars(16); // MU mf(5) = pars(17); // MS mf(2) = pars(18); // MC + mcMc = pars(19); + mu2GeV = pars(20); + md2GeV = pars(21); + ms2GeV = pars(22); } Eigen::ArrayXd QedQcd::display_input() const { - Eigen::ArrayXd pars(19); + Eigen::ArrayXd pars(23); pars(0) = a(1); pars(1) = gfermi; @@ -731,6 +764,10 @@ Eigen::ArrayXd QedQcd::display_input() const pars(16) = mf(1); // MU pars(17) = mf(5); // MS pars(18) = mf(2); // MC + pars(19) = mcMc; + pars(20) = mu2GeV; + pars(21) = md2GeV; + pars(22) = ms2GeV; return pars; } diff --git a/src/lowe.h b/src/lowe.h index 3d0edf8cf..0e3003aec 100644 --- a/src/lowe.h +++ b/src/lowe.h @@ -53,6 +53,7 @@ typedef enum {ALPHA=1, ALPHAS} leGauge; enum QedQcd_input_parmeters : unsigned { alpha_em_MSbar_at_MZ, GFermi, alpha_s_MSbar_at_MZ, MZ_pole, mb_mb, MT_pole, MTau_pole, MMuon_pole, MElectron_pole, Mv3_pole, MW_pole, ME_pole, Mv1_pole, MM_pole, Mv2_pole, MD, MU, MS, MC, + mc_mc, Mu2GeVInput, Md2GeVInput, Ms2GeVInput, NUMBER_OF_LOW_ENERGY_INPUT_PARAMETERS}; extern const char* QedQcd_input_parmeter_names[NUMBER_OF_LOW_ENERGY_INPUT_PARAMETERS]; @@ -77,6 +78,8 @@ class QedQcd: public RGE double mwPole; ///< W boson pole mass double mzPole; ///< Z boson pole mass double gfermi; ///< Fermi constant + double mcMc; ///< mc(mc) + double mu2GeV, md2GeV, ms2GeV; flexiblesusy::CKM_parameters ckm; ///< CKM parameters (in the MS-bar scheme at MZ) flexiblesusy::PMNS_parameters pmns; ///< PMNS parameters (in the MS-bar scheme at MZ) @@ -94,6 +97,10 @@ class QedQcd: public RGE void setPoleMmuon(double m) { mmuonPole = m; } ///< set pole muon mass void setPoleMel(double m) { melPole = m; } ///< set pole electron mass void setMbMb(double mb) { mbMb = mb; }; ///< set mb(mb) + void setMcMc(double mc) { mcMc = mc; } ///< set mc(mc) + void setMu2GeV(double mu) { mu2GeV = mu; } ///< set mu(2 GeV) + void setMd2GeV(double md) { md2GeV = md; } ///< set md(2 GeV) + void setMs2GeV(double ms) { ms2GeV = ms; } ///< set ms(2 GeV) void setPoleMW(double mw) { mwPole = mw; } ///< set W boson pole mass void setPoleMZ(double mz) { mzPole = mz; } ///< set Z boson pole mass /// sets a running quark mass @@ -145,6 +152,14 @@ class QedQcd: public RGE static std::vector display_input_parameter_names(); /// Returns mb(mb) MSbar double displayMbMb() const { return mbMb; } + /// Returns mc(mc) MSbar + double displayMcMc() const { return mcMc; } + /// Returns mu(2 GeV) + double displayMu2GeV() const { return mu2GeV; } + /// Returns md(2 GeV) + double displayMd2GeV() const { return md2GeV; } + /// Returns ms(2 GeV) + double displayMs2GeV() const { return ms2GeV; } /// returns CKM parameters flexiblesusy::CKM_parameters displayCKM() const { return ckm; } /// Returns real CKM matrix @@ -216,6 +231,7 @@ inline QedQcd::QedQcd(const QedQcd &m) : RGE(), a(m.a), mf(m.mf), mnu(m.mnu), mtPole(m.mtPole), mbPole(m.mbPole), mbMb(m.mbMb), mtauPole(m.mtauPole), mmuonPole(m.mmuonPole), melPole(m.melPole), mwPole(m.mwPole), mzPole(m.mzPole), gfermi(m.gfermi), ckm(m.ckm), pmns(m.pmns) + , mcMc(m.mcMc), mu2GeV(m.mu2GeV), md2GeV(m.md2GeV), ms2GeV(m.ms2GeV) { setPars(11); setMu(m.displayMu()); diff --git a/src/slha_io.cpp b/src/slha_io.cpp index 95ed8e318..88a4b3fc1 100644 --- a/src/slha_io.cpp +++ b/src/slha_io.cpp @@ -520,6 +520,7 @@ void SLHA_io::process_sminputs_tuple(softsusy::QedQcd& qedqcd, int key, double v break; case 4: qedqcd.setPoleMZ(value); + qedqcd.setMu(value); softsusy::MZ = value; break; case 5: @@ -555,15 +556,19 @@ void SLHA_io::process_sminputs_tuple(softsusy::QedQcd& qedqcd, int key, double v break; case 21: qedqcd.setMass(mDown, value); + qedqcd.setMd2GeV(value); break; case 22: qedqcd.setMass(mUp, value); + qedqcd.setMu2GeV(value); break; case 23: qedqcd.setMass(mStrange, value); + qedqcd.setMs2GeV(value); break; case 24: qedqcd.setMass(mCharm, value); + qedqcd.setMcMc(value); break; default: WARNING("Unrecognized entry in block SMINPUTS: " << key); diff --git a/test/test_lowe.cpp b/test/test_lowe.cpp index 67dfece9d..c5f18f191 100644 --- a/test/test_lowe.cpp +++ b/test/test_lowe.cpp @@ -35,3 +35,24 @@ BOOST_AUTO_TEST_CASE( test_to ) BOOST_MESSAGE(lowe_Mz); BOOST_MESSAGE(lowe_Mz_new); } + +BOOST_AUTO_TEST_CASE( test_to_recall ) +{ + QedQcd lowe_Mz; + lowe_Mz.setPoleMt(173.5); + lowe_Mz.setAlpha(ALPHAS, 0.118); + lowe_Mz.setMu(lowe_Mz.displayPoleMZ()); + lowe_Mz.to(100.); + + QedQcd lowe_Mz_new(lowe_Mz); + lowe_Mz.setMu(100.); + lowe_Mz_new.to(100.); + + BOOST_CHECK_LT(flexiblesusy::MaxRelDiff(lowe_Mz.display_input(), lowe_Mz_new.display_input()), 5e-4); + + BOOST_MESSAGE(lowe_Mz.display_input().transpose()); + BOOST_MESSAGE(lowe_Mz_new.display_input().transpose()); + + BOOST_MESSAGE(lowe_Mz); + BOOST_MESSAGE(lowe_Mz_new); +}