standard_model_two_scale_high_scale_spectrum_generator.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_spectrum_generator.hpp"
#include "@ModelName@_two_scale_convergence_tester.hpp"
#include "@ModelName@_two_scale_initial_guesser.hpp"
#include "standard_model_two_scale_convergence_tester.hpp"
#include "two_scale_composite_convergence_tester.hpp"

#include "error.hpp"
#include "lowe.h"
#include "numerics2.hpp"
#include "two_scale_running_precision.hpp"

#include <algorithm>
#include <limits>

namespace flexiblesusy {

double @ModelName@_spectrum_generator<Two_scale>::get_pole_mass_scale() const
{
   return
      this->settings.get(Spectrum_generator_settings::pole_mass_scale) != 0. ?
      this->settings.get(Spectrum_generator_settings::pole_mass_scale) : get_susy_scale();
}

double @ModelName@_spectrum_generator<Two_scale>::get_eft_pole_mass_scale() const
{
   double Q_higgs = this->settings.get(Spectrum_generator_settings::eft_pole_mass_scale);

   if (Q_higgs == 0.) {
      const double Mt = low_scale_constraint.get_sm_parameters().displayPoleMt();
      Q_higgs = std::min(get_susy_scale(), Mt);
   }

   return Q_higgs;
}

/**
 * @brief Run's the RG solver with the given input parameters
 *
 * This function sets up the RG solver using a high-scale, susy-scale
 * and low-scale constraint.  Afterwards the solver is run until
 * convergence is reached or an error occours.  Finally the particle
 * spectrum (pole masses) is calculated.
 *
 * @param qedqcd Standard Model input parameters
 * @param input model input parameters
 */
void @ModelName@_spectrum_generator<Two_scale>::run_except(const softsusy::QedQcd& qedqcd,
                                const @ModelName@_input_parameters& input)
{
   @ModelName@<Two_scale>& model = this->model;
   model.clear();
   model.set_input_parameters(input);
   model.do_calculate_sm_pole_masses(this->settings.get(Spectrum_generator_settings::calculate_sm_masses));
   model.do_calculate_bsm_pole_masses(this->settings.get(Spectrum_generator_settings::calculate_bsm_masses));
   model.do_force_output(this->settings.get(Spectrum_generator_settings::force_output));
   model.set_loops(this->settings.get(Spectrum_generator_settings::beta_loop_order));
   model.set_thresholds(this->settings.get(Spectrum_generator_settings::threshold_corrections_loop_order));
   model.set_zero_threshold(this->settings.get(Spectrum_generator_settings::beta_zero_threshold));

   eft.clear();
   eft.do_force_output(this->settings.get(Spectrum_generator_settings::force_output));
   eft.set_loops(this->settings.get(Spectrum_generator_settings::beta_loop_order));
   eft.set_thresholds(this->settings.get(Spectrum_generator_settings::threshold_corrections_loop_order));
   eft.set_zero_threshold(this->settings.get(Spectrum_generator_settings::beta_zero_threshold));
   eft.set_pole_mass_loop_order(this->model.get_pole_mass_loop_order());
   eft.set_ewsb_loop_order(this->model.get_ewsb_loop_order());
   eft.set_ewsb_iteration_precision(this->model.get_ewsb_iteration_precision());
   eft.set_number_of_ewsb_iterations(this->model.get_number_of_ewsb_iterations());
   eft.set_number_of_mass_iterations(this->model.get_number_of_mass_iterations());
   eft.set_two_loop_corrections(this->model.get_two_loop_corrections());

   high_scale_constraint.clear();
   susy_scale_constraint.clear();
   low_scale_constraint .clear();

   // needed for constraint::initialize()
   high_scale_constraint.set_model(&model);
   susy_scale_constraint.set_model(&model);
   low_scale_constraint .set_model(&eft);

   low_scale_constraint .set_sm_parameters(qedqcd);

   const unsigned index = this->settings.get(Spectrum_generator_settings::eft_higgs_index);
   const auto scale_getter = [this] () { return susy_scale_constraint.get_scale(); };

   matching_up.set_models(&eft, &model);
   matching_up.set_scale(scale_getter);
   matching_up.set_scale(this->settings.get(Spectrum_generator_settings::eft_matching_scale));
   matching_up.set_loop_order(this->settings.get(Spectrum_generator_settings::eft_matching_loop_order_up));
   matching_up.set_higgs_index(index);

   matching_down.set_models(&model, &eft);
   matching_down.set_scale(scale_getter);
   matching_down.set_scale(this->settings.get(Spectrum_generator_settings::eft_matching_scale));
   matching_down.set_loop_order(this->settings.get(Spectrum_generator_settings::eft_matching_loop_order_down));
   matching_down.set_higgs_index(index);

   high_scale_constraint.initialize();
   susy_scale_constraint.initialize();
   low_scale_constraint .initialize();

   // convergence tester for @ModelName@
   @ModelName@_convergence_tester<Two_scale> model_ct(
      &model, this->settings.get(Spectrum_generator_settings::precision),
      [this](){ return get_pole_mass_scale(); });

   // convergence tester for Standard Model
   standard_model::Standard_model_convergence_tester<Two_scale> eft_ct(
      &eft, this->settings.get(Spectrum_generator_settings::precision),
      [this](){ return get_eft_pole_mass_scale(); });

   if (this->settings.get(Spectrum_generator_settings::max_iterations) > 0) {
      model_ct.set_max_iterations(
         this->settings.get(Spectrum_generator_settings::max_iterations));
      eft_ct.set_max_iterations(
         this->settings.get(Spectrum_generator_settings::max_iterations));
   }

   Composite_convergence_tester<Two_scale> cct;
   cct.add_convergence_tester(&model_ct);
   cct.add_convergence_tester(&eft_ct);

   @ModelName@_standard_model_initial_guesser<Two_scale> initial_guesser(&model, &eft, qedqcd,
                                                                         low_scale_constraint,
                                                                         susy_scale_constraint,
                                                                         high_scale_constraint);

   Two_scale_increasing_precision precision(
      10.0, this->settings.get(Spectrum_generator_settings::precision));

   solver.reset();
   solver.set_convergence_tester(&cct);
   solver.set_running_precision(&precision);
   solver.set_initial_guesser(&initial_guesser);
   solver.add(&low_scale_constraint, &eft);
   solver.add(&matching_up, &eft, &model);
   solver.add(&susy_scale_constraint, &model);
   solver.add(&high_scale_constraint, &model);
   solver.add(&susy_scale_constraint, &model);
   solver.add(&matching_down, &model, &eft);

   this->reached_precision = std::numeric_limits<double>::infinity();

   solver.solve();

   // impose low-scale constraint one last time
   eft.run_to(low_scale_constraint.get_scale());
   low_scale_constraint.apply();

   this->reached_precision = std::max(model_ct.get_current_accuracy(),
                                      eft_ct.get_current_accuracy());

   calculate_spectrum();

   // copy calculated W pole mass
   model.get_physical().M@VectorW@
      = low_scale_constraint.get_sm_parameters().displayPoleMW();

   // run to output scale (if scale > 0)
   if (!is_zero(this->parameter_output_scale)) {
      model.run_to(this->parameter_output_scale);
      eft.run_to(this->parameter_output_scale);
   }
}

/**
 * Create a text file which contains the values of all model
 * parameters at all scales between the low-scale and the high-scale.
 *
 * @param filename name of output file
 */
void @ModelName@_spectrum_generator<Two_scale>::write_running_couplings(
   const std::string& filename) const
{
   @ModelName@_spectrum_generator_interface<Two_scale>::write_running_couplings(
      filename, get_low_scale(), get_high_scale());
}

void @ModelName@_spectrum_generator<Two_scale>::calculate_spectrum()
{
   this->model.run_to(get_pole_mass_scale());
   this->model.solve_ewsb();
   this->model.calculate_spectrum();

   eft.run_to(get_eft_pole_mass_scale());

   // computation of pole mass spectrum in the SM
   const unsigned eft_pole_loops = eft.get_pole_mass_loop_order();
   const unsigned eft_ewsb_loops = eft.get_ewsb_loop_order();

   if (eft_pole_loops > 1) {
      WARNING("Pole mass loop order " << eft_pole_loops
              << " in the EFT is currently not supported!  I'm using 1-loop.");
      eft.set_pole_mass_loop_order(1);
   }
   if (eft_ewsb_loops > 1) {
      WARNING("EWSB loop order " << eft_ewsb_loops
              << " in the EFT is currently not supported!  I'm using 1-loop.");
      eft.set_ewsb_loop_order(1);
   }

   eft.calculate_DRbar_masses();
   eft.solve_ewsb();
   eft.calculate_spectrum();

   eft.set_pole_mass_loop_order(eft_pole_loops);
   eft.set_ewsb_loop_order(eft_ewsb_loops);

   const unsigned index = this->settings.get(Spectrum_generator_settings::eft_higgs_index);

   this->model.get_physical().M@HiggsBoson_index@ = eft.get_physical().Mhh;
   this->model.get_physical().M@VectorZ_0@ = eft.get_physical().MVZ;
   this->model.get_physical().M@VectorW_0@ = eft.get_physical().MVWp;
@fillSMFermionPoleMasses@
   if (eft.get_problems().is_running_tachyon(standard_model_info::hh))
      this->model.get_problems().flag_running_tachyon(@ModelName@_info::@HiggsBoson@);
   if (eft.get_problems().is_pole_tachyon(standard_model_info::hh))
      this->model.get_problems().flag_pole_tachyon(@ModelName@_info::@HiggsBoson@);
   if (eft.get_problems().is_running_tachyon(standard_model_info::VZ))
      this->model.get_problems().flag_running_tachyon(@ModelName@_info::@VectorZ@);
   if (eft.get_problems().is_pole_tachyon(standard_model_info::VZ))
      this->model.get_problems().flag_pole_tachyon(@ModelName@_info::@VectorZ@);
   if (eft.get_problems().is_running_tachyon(standard_model_info::VWp))
      this->model.get_problems().flag_running_tachyon(@ModelName@_info::@VectorW@);
   if (eft.get_problems().is_pole_tachyon(standard_model_info::VWp))
      this->model.get_problems().flag_pole_tachyon(@ModelName@_info::@VectorW@);
}

} // namespace flexiblesusy