two_scale_low_scale_constraint.cpp.in

// ====================================================================
// 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/>.
// ====================================================================

// File generated at @DateAndTime@

#include "@ModelName@_two_scale_low_scale_constraint.hpp"
#include "@ModelName@_two_scale_model.hpp"
#include "@ModelName@_info.hpp"
#include "wrappers.hpp"
#include "logger.hpp"
#include "error.hpp"
#include "ew_input.hpp"
#include "gsl_utils.hpp"
#include "minimizer.hpp"
#include "raii.hpp"
#include "root_finder.hpp"
#include "threshold_loop_functions.hpp"
#include "weinberg_angle.hpp"

#include <cmath>
#include <limits>

namespace flexiblesusy {

#define DERIVEDPARAMETER(p) model->p()
#define INPUTPARAMETER(p) model->get_input().p
#define MODELPARAMETER(p) model->get_##p()
#define PHASE(p) model->get_##p()
#define BETAPARAMETER(p) beta_functions.get_##p()
#define BETAPARAMETER1(l,p) beta_functions_##l##L.get_##p()
#define BETA(p) beta_##p
#define BETA1(l,p) beta_##l##L_##p
#define LowEnergyGaugeCoupling(i) new_g##i
#define LowEnergyConstant(p) Electroweak_constants::p
#define MZPole qedqcd.displayPoleMZ()
#define STANDARDDEVIATION(p) Electroweak_constants::Error_##p
#define Pole(p) model->get_physical().p
#define SCALE model->get_scale()
#define MODEL model
#define MODELCLASSNAME @ModelName@<Two_scale>
#define CKM ckm
#define PMNS pmns
#define THETAW theta_w
#define THRESHOLD static_cast<int>(model->get_thresholds())
#define ALPHA_EM_DRBAR alpha_em_drbar
#define CALCULATE_DRBAR_MASSES() model->calculate_DRbar_masses()

@ModelName@_low_scale_constraint<Two_scale>::@ModelName@_low_scale_constraint(
   @ModelName@<Two_scale>* model_, const softsusy::QedQcd& qedqcd_)
   : model(model_)
   , qedqcd(qedqcd_)
{
   initialize();
}

@ModelName@_low_scale_constraint<Two_scale>::~@ModelName@_low_scale_constraint()
{
}

void @ModelName@_low_scale_constraint<Two_scale>::apply()
{
   check_model_ptr();

@temporarySetting@

   model->calculate_DRbar_masses();
   update_scale();
   qedqcd.run_to(scale, 1.0e-5);
   calculate_DRbar_gauge_couplings();

@applyConstraint@

   if (model->get_thresholds())
      qedqcd.setPoleMW(recalculate_mw_pole(qedqcd.displayPoleMW()));
}

const Eigen::Matrix<std::complex<double>,3,3>& @ModelName@_low_scale_constraint<Two_scale>::get_ckm()
{
   return ckm;
}

const Eigen::Matrix<std::complex<double>,3,3>& @ModelName@_low_scale_constraint<Two_scale>::get_pmns()
{
   return pmns;
}

double @ModelName@_low_scale_constraint<Two_scale>::get_scale() const
{
   return scale;
}

double @ModelName@_low_scale_constraint<Two_scale>::get_initial_scale_guess() const
{
   return initial_scale_guess;
}

void @ModelName@_low_scale_constraint<Two_scale>::set_model(Model* model_)
{
   model = cast_model<@ModelName@<Two_scale>*>(model_);
}

void @ModelName@_low_scale_constraint<Two_scale>::set_sm_parameters(
   const softsusy::QedQcd& qedqcd_)
{
   qedqcd = qedqcd_;
}

const softsusy::QedQcd& @ModelName@_low_scale_constraint<Two_scale>::get_sm_parameters() const
{
   return qedqcd;
}

void @ModelName@_low_scale_constraint<Two_scale>::clear()
{
   scale = 0.;
   initial_scale_guess = 0.;
   model = nullptr;
   qedqcd = softsusy::QedQcd();
   ckm.setIdentity();
   pmns.setIdentity();
   neutrinoDRbar.setZero();
   MWDRbar = 0.;
   MZDRbar = 0.;
   AlphaS = 0.;
   EDRbar = 0.;
   ThetaWDRbar = 0.;
   new_g1 = 0.;
   new_g2 = 0.;
   new_g3 = 0.;
   self_energy_w_at_mw = 0.;
}

void @ModelName@_low_scale_constraint<Two_scale>::initialize()
{
   check_model_ptr();

@scaleGuess@
   scale = initial_scale_guess;

   ckm = qedqcd.get_complex_ckm();
   pmns = qedqcd.get_complex_pmns();
   neutrinoDRbar = Eigen::Matrix<double,3,3>::Zero();
   MWDRbar = 0.;
   MZDRbar = 0.;
   AlphaS = 0.;
   EDRbar = 0.;
   ThetaWDRbar = 0.;
   new_g1 = 0.;
   new_g2 = 0.;
   new_g3 = 0.;
   self_energy_w_at_mw = 0.;
}

void @ModelName@_low_scale_constraint<Two_scale>::update_scale()
{
   check_model_ptr();

@calculateScale@
@restrictScale@
}

void @ModelName@_low_scale_constraint<Two_scale>::calculate_threshold_corrections()
{
   check_model_ptr();

   if (qedqcd.get_scale() != get_scale())
      throw SetupError("Error: low-energy data"
          " set is not defined at the same scale as the low-energy"
          " constraint.  You need to run the low-energy data set to this"
          " scale!");

   const double alpha_em = qedqcd.displayAlpha(softsusy::ALPHA);
   const double alpha_s  = qedqcd.displayAlpha(softsusy::ALPHAS);
   const double mw_pole  = qedqcd.displayPoleMW();
   const double mz_pole  = qedqcd.displayPoleMZ();

   double delta_alpha_em = 0.;
   double delta_alpha_s  = 0.;

   if (model->get_thresholds()) {
      delta_alpha_em = calculate_delta_alpha_em(alpha_em);
      delta_alpha_s  = calculate_delta_alpha_s(alpha_s);
   }

   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);

   // interface variables
   MZDRbar = mz_pole;
   MWDRbar = mw_pole;

   if (model->get_thresholds()) {
      MZDRbar = model->calculate_M@VectorZ@_DRbar(mz_pole);
      MWDRbar = model->calculate_M@VectorW@_DRbar(mw_pole);
   }

   AlphaS = alpha_s_drbar;
   EDRbar = e_drbar;
   ThetaWDRbar = calculate_theta_w(alpha_em_drbar);
}

double @ModelName@_low_scale_constraint<Two_scale>::calculate_theta_w(double alpha_em_drbar)
{
   check_model_ptr();

   double theta_w = 0.;

@calculateThetaW@

   return theta_w;
}

void @ModelName@_low_scale_constraint<Two_scale>::calculate_DRbar_gauge_couplings()
{
   check_model_ptr();
   calculate_threshold_corrections();

@calculateGaugeCouplings@

   if (IsFinite(new_g1)) {
      model->get_problems().unflag_non_perturbative_parameter(@ModelName@_info::@hyperchargeCoupling@);
   } else {
      model->get_problems().flag_non_perturbative_parameter(
         @ModelName@_info::@hyperchargeCoupling@, new_g1, get_scale());
      new_g1 = Electroweak_constants::g1;
   }

   if (IsFinite(new_g2)) {
      model->get_problems().unflag_non_perturbative_parameter(@ModelName@_info::@leftCoupling@);
   } else {
      model->get_problems().flag_non_perturbative_parameter(
         @ModelName@_info::@leftCoupling@, new_g2, get_scale());
      new_g2 = Electroweak_constants::g2;
   }
}

double @ModelName@_low_scale_constraint<Two_scale>::calculate_delta_alpha_em(double alphaEm) const
{
   check_model_ptr();

   const double currentScale = model->get_scale();
@calculateDeltaAlphaEm@
}

double @ModelName@_low_scale_constraint<Two_scale>::calculate_delta_alpha_s(double alphaS) const
{
   check_model_ptr();

   const double currentScale = model->get_scale();
@calculateDeltaAlphaS@
}

void @ModelName@_low_scale_constraint<Two_scale>::calculate_DRbar_yukawa_couplings()
{
   calculate_@UpYukawa@_DRbar();
   calculate_@DownYukawa@_DRbar();
   calculate_@ElectronYukawa@_DRbar();
}

void @ModelName@_low_scale_constraint<Two_scale>::calculate_@UpYukawa@_DRbar()
{
   check_model_ptr();

   Eigen::Matrix<std::complex<double>,3,3> upQuarksDRbar(ZEROMATRIXCOMPLEX(3,3));
   upQuarksDRbar(0,0)      = qedqcd.displayMass(softsusy::mUp);
   upQuarksDRbar(1,1)      = qedqcd.displayMass(softsusy::mCharm);
   upQuarksDRbar(2,2)      = qedqcd.displayPoleMt();

   if (model->get_thresholds()) {
@calculateDRbarMassTop@
   }

@setDRbarUpQuarkYukawaCouplings@
}

void @ModelName@_low_scale_constraint<Two_scale>::calculate_@DownYukawa@_DRbar()
{
   check_model_ptr();

   Eigen::Matrix<std::complex<double>,3,3> downQuarksDRbar(ZEROMATRIXCOMPLEX(3,3));
   downQuarksDRbar(0,0)   = qedqcd.displayMass(softsusy::mDown);
   downQuarksDRbar(1,1)   = qedqcd.displayMass(softsusy::mStrange);
   downQuarksDRbar(2,2)   = qedqcd.displayMass(softsusy::mBottom);

   if (model->get_thresholds()) {
@calculateDRbarMassBottom@
   }

@setDRbarDownQuarkYukawaCouplings@
}

void @ModelName@_low_scale_constraint<Two_scale>::calculate_@ElectronYukawa@_DRbar()
{
   check_model_ptr();

   Eigen::Matrix<std::complex<double>,3,3> downLeptonsDRbar(ZEROMATRIXCOMPLEX(3,3));
   downLeptonsDRbar(0,0) = qedqcd.displayPoleMel();
   downLeptonsDRbar(1,1) = qedqcd.displayPoleMmuon();
   downLeptonsDRbar(2,2) = qedqcd.displayPoleMtau();

   if (model->get_thresholds()) {
@calculateDRbarMassElectron@
@calculateDRbarMassMuon@
@calculateDRbarMassTau@
   }

@setDRbarElectronYukawaCouplings@
}

void @ModelName@_low_scale_constraint<Two_scale>::calculate_MNeutrino_DRbar()
{
   neutrinoDRbar.setZero();
   neutrinoDRbar(0,0) = qedqcd.displayNeutrinoPoleMass(1);
   neutrinoDRbar(1,1) = qedqcd.displayNeutrinoPoleMass(2);
   neutrinoDRbar(2,2) = qedqcd.displayNeutrinoPoleMass(3);
}

/**
 * Recalculates the W boson pole mass using the new gauge couplings.
 *
 * @param mw_pole current W pole mass
 * @return new W pole mass
 */
double @ModelName@_low_scale_constraint<Two_scale>::recalculate_mw_pole(double mw_pole)
{
   check_model_ptr();

@recalculateMWPole@
}

void @ModelName@_low_scale_constraint<Two_scale>::check_model_ptr() const
{
   if (!model)
      throw SetupError("@ModelName@_low_scale_constraint<Two_scale>: "
                       "model pointer is zero!");
}

} // namespace flexiblesusy