standard_model_high_scale_spectrum_generator.hpp.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@

#ifndef @ModelName@_STANDARD_MODEL_SPECTRUM_GENERATOR_H
#define @ModelName@_STANDARD_MODEL_SPECTRUM_GENERATOR_H

#include "@ModelName@_spectrum_generator_interface.hpp"
#include "@ModelName@_two_scale_high_scale_constraint.hpp"
#include "@ModelName@_two_scale_susy_scale_constraint.hpp"
#include "@ModelName@_two_scale_convergence_tester.hpp"
#include "@ModelName@_two_scale_initial_guesser.hpp"
#include "@ModelName@_standard_model_matching.hpp"
#include "@ModelName@_input_parameters.hpp"
#include "@ModelName@_standard_model_two_scale_matching.hpp"
#include "standard_model_two_scale_model.hpp"
#include "standard_model_two_scale_low_scale_constraint.hpp"
#include "standard_model_two_scale_convergence_tester.hpp"
#include "two_scale_composite_convergence_tester.hpp"

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

#include <algorithm>
#include <limits>

namespace flexiblesusy {

template <class T>
class @ModelName@_spectrum_generator
   : public @ModelName@_spectrum_generator_interface<T> {
public:
   @ModelName@_spectrum_generator()
      : @ModelName@_spectrum_generator_interface<T>()
      , eft()
      , high_scale(0.)
      , susy_scale(0.)
      , low_scale(0.)
   {}
   virtual ~@ModelName@_spectrum_generator() {}

   const standard_model::StandardModel<T>& get_eft() const { return eft; }
   standard_model::StandardModel<T>& get_eft() { return eft; }
   double get_high_scale() const { return high_scale; }
   double get_susy_scale() const { return susy_scale; }
   double get_low_scale()  const { return low_scale;  }

   virtual void run(const softsusy::QedQcd&, const @ModelName@_input_parameters&);
   void write_running_couplings(const std::string& filename = "@ModelName@_rgflow.dat") const;

private:
   standard_model::StandardModel<T> eft;
   double high_scale, susy_scale, low_scale;
};

/**
 * @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
 */
template <class T>
void @ModelName@_spectrum_generator<T>::run(const softsusy::QedQcd& qedqcd,
                                const @ModelName@_input_parameters& input)
{
   @ModelName@<T>& 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());

   @ModelName@_high_scale_constraint<T> high_scale_constraint;
   @ModelName@_susy_scale_constraint<T> susy_scale_constraint;
   standard_model::Standard_model_low_scale_constraint<T>  low_scale_constraint;

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

   @ModelName@_standard_model_Matching<T> matching;
   matching.set_models(&eft, &model);
   matching.set_constraint(&susy_scale_constraint);
   matching.set_scale(this->settings.get(Spectrum_generator_settings::eft_matching_scale));
   matching.set_loop_order_up(this->settings.get(Spectrum_generator_settings::eft_matching_loop_order_up));
   matching.set_loop_order_down(this->settings.get(Spectrum_generator_settings::eft_matching_loop_order_down));
   matching.set_higgs_index(index);

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

   std::vector<Constraint<T>*> model_constraints(2);
   model_constraints[0] = &susy_scale_constraint;
   model_constraints[1] = &high_scale_constraint;

   std::vector<Constraint<T>*> eft_constraints(1);
   eft_constraints[0] = &low_scale_constraint;

   // convergence tester for @ModelName@
   @ModelName@_convergence_tester<T> model_ct(
      &model, this->settings.get(Spectrum_generator_settings::precision));

   // convergence tester for Standard Model
   standard_model::Standard_model_convergence_tester<T> eft_ct(
      &eft, this->settings.get(Spectrum_generator_settings::precision));

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

   RGFlow<T> solver;
   solver.set_convergence_tester(&cct);
   solver.set_running_precision(&precision);
   solver.set_initial_guesser(&initial_guesser);
   solver.add_model(&eft, &matching, eft_constraints);
   solver.add_model(&model, model_constraints);

   high_scale = susy_scale = low_scale = 0.;
   this->reached_precision = std::numeric_limits<double>::infinity();

   try {
      solver.solve();
      high_scale = high_scale_constraint.get_scale();
      susy_scale = susy_scale_constraint.get_scale();
      low_scale  = low_scale_constraint.get_scale();
      this->reached_precision = std::max(model_ct.get_current_accuracy(),
                                         eft_ct.get_current_accuracy());

      const double mass_scale =
         this->settings.get(Spectrum_generator_settings::pole_mass_scale) != 0. ?
         this->settings.get(Spectrum_generator_settings::pole_mass_scale) : susy_scale;

      model.run_to(mass_scale);
      model.solve_ewsb();
      model.calculate_spectrum();

      // start logarithmic resummation
      {
         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(susy_scale, Mt);
         }

         eft.run_to(Q_higgs);

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

         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_tachyon(standard_model_info::hh))
            this->model.get_problems().flag_tachyon(@ModelName@_info::@HiggsBoson@);
         if (eft.get_problems().is_tachyon(standard_model_info::VZ))
            this->model.get_problems().flag_tachyon(@ModelName@_info::@VectorZ@);
         if (eft.get_problems().is_tachyon(standard_model_info::VWp))
            this->model.get_problems().flag_tachyon(@ModelName@_info::@VectorW@);
      }

      // 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);
      }
   } catch (...) {
      this->translate_exception_to_problem(model);
   }
}

/**
 * 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
 */
template <class T>
void @ModelName@_spectrum_generator<T>::write_running_couplings(
   const std::string& filename) const
{
   @ModelName@_spectrum_generator_interface<T>::write_running_couplings(filename, low_scale, high_scale);
}

} // namespace flexiblesusy

#endif