Skip to content

Commit

Permalink
make to() re-callable
Browse files Browse the repository at this point in the history 
by adding fixed input parameters
  • Loading branch information
Alexander Voigt authored and Alexander Voigt committed Apr 20, 2016
1 parent 7c2ba80 commit 3a585aa
Show file tree
Hide file tree
Showing 4 changed files with 100 additions and 21 deletions.
79 changes: 58 additions & 21 deletions src/lowe.cpp
View file
Expand Up @@ -9,15 +9,14 @@
#include "lowe.h"
#include "ew_input.hpp"
#include "wrappers.hpp"
#include "error.hpp"

namespace softsusy {

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!
Expand All @@ -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'
Expand Down Expand Up @@ -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());
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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;
}
Expand Down
16 changes: 16 additions & 0 deletions src/lowe.h
View file
Expand Up @@ -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];
Expand All @@ -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)

Expand All @@ -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
Expand Down Expand Up @@ -145,6 +152,14 @@ class QedQcd: public RGE
static std::vector<std::string> 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
Expand Down Expand Up @@ -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());
Expand Down
5 changes: 5 additions & 0 deletions src/slha_io.cpp
View file
Expand Up @@ -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:
Expand Down Expand Up @@ -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);
Expand Down
21 changes: 21 additions & 0 deletions test/test_lowe.cpp
View file
Expand Up @@ -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);
}

0 comments on commit 3a585aa

Please sign in to comment.