ClassPotentialOrigin.h

Go to the documentation of this file.

// Copyright (C) 2018  Philipp Basler and Margarete Mühlleitner
// SPDX-FileCopyrightText: 2021 Philipp Basler, Margarete Mühlleitner and Jonas
// Müller
//
// SPDX-License-Identifier: GPL-3.0-or-later
 
#ifndef CLASSPOTENTIALORIGIN_H_
#define CLASSPOTENTIALORIGIN_H_
 
#include "Eigen/Eigenvalues"
#include <BSMPT/minimizer/Minimizer.h>
#include <BSMPT/models/IncludeAllModels.h>
#include <BSMPT/models/SMparam.h>
#include <BSMPT/utility/settings.h>
#include <iostream>
#include <vector>
 
namespace BSMPT
{
 
const double C_PT = 0;
 
const double C_threshold = 1e-4;
static constexpr double C_CWcbFermion = 1.5;
const double C_CWcbGB = 5.0 / 6.0;
const double C_CWcbHiggs = 1.5;
 
enum class Order
{
  TreeLevel,
  OneLoop
};
 
class Class_Potential_Origin
{
public:
  const ISMConstants SMConstants;
 
protected:
  bool UseTreeLevel = false;
 
  double scale;
 
  std::size_t nPar = 0;
  std::size_t nParCT = 0;
 
  std::vector<double> parStored;
 
  std::vector<double> parCTStored;
 
  ModelID::ModelIDs Model = ModelID::ModelIDs::NotSet;
  std::size_t NNeutralHiggs = 0;
  std::size_t NChargedHiggs = 0;
  std::size_t NHiggs = NNeutralHiggs + NChargedHiggs;
  std::size_t NGauge = 4;
  std::size_t NQuarks = 12;
 
  std::size_t NColour = 3;
 
  std::size_t NLepton = 9;
 
  std::size_t nVEV = 0;
 
  std::vector<double> VEVSymmetric;
  std::vector<double> vevTree;
  std::vector<double> vevTreeMin;
  std::vector<double> TripleHiggsCorrectionsCW;
  std::vector<std::vector<std::vector<double>>>
      TripleHiggsCorrectionsCWPhysical;
  std::vector<std::vector<std::vector<double>>>
      TripleHiggsCorrectionsTreePhysical;
  std::vector<std::vector<std::vector<double>>>
      TripleHiggsCorrectionsCTPhysical;
 
  bool SetCurvatureDone = false;
  bool CalcCouplingsDone = false;
  bool CalculatedTripleCopulings = false;
 
  int InputLineNumber = -1;
 
  bool UseIndexCol = false;
 
  std::vector<double> HiggsVev;
  std::vector<double> Curvature_Higgs_L1;
  std::vector<std::vector<double>> Curvature_Higgs_L2;
  std::vector<std::vector<std::vector<double>>> Curvature_Higgs_L3;
  std::vector<std::vector<std::vector<std::vector<double>>>> Curvature_Higgs_L4;
  std::vector<double> Curvature_Higgs_CT_L1;
  std::vector<std::vector<double>> Curvature_Higgs_CT_L2;
  std::vector<std::vector<std::vector<double>>> Curvature_Higgs_CT_L3;
  std::vector<std::vector<std::vector<std::vector<double>>>>
      Curvature_Higgs_CT_L4;
  std::vector<std::vector<std::vector<std::vector<double>>>>
      Curvature_Gauge_G2H2;
  std::vector<std::vector<std::vector<std::complex<double>>>>
      Curvature_Quark_F2H1;
  std::vector<std::vector<std::complex<double>>> Curvature_Quark_F2;
  std::vector<std::vector<std::vector<std::complex<double>>>>
      Curvature_Lepton_F2H1;
  std::vector<std::vector<std::complex<double>>> Curvature_Lepton_F2;
  std::vector<double> MassSquaredHiggs;
  std::vector<double> MassSquaredGauge;
  std::vector<double> MassSquaredQuark;
  std::vector<double> MassSquaredLepton;
  std::vector<std::vector<double>> HiggsRotationMatrix;
  std::vector<std::vector<double>> HiggsRotationMatrixEnsuredConvention;
  std::vector<std::vector<std::vector<std::vector<double>>>>
      Couplings_Higgs_Quartic;
  std::vector<std::vector<std::vector<double>>> Couplings_Higgs_Triple;
  std::vector<std::vector<std::vector<std::vector<double>>>>
      Couplings_Gauge_Higgs_22;
  std::vector<std::vector<std::vector<double>>> Couplings_Gauge_Higgs_21;
  std::vector<std::vector<std::vector<std::vector<std::complex<double>>>>>
      Couplings_Quark_Higgs_22;
  std::vector<std::vector<std::vector<std::complex<double>>>>
      Couplings_Quark_Higgs_21;
  std::vector<std::vector<std::vector<std::vector<std::complex<double>>>>>
      Couplings_Lepton_Higgs_22;
  std::vector<std::vector<std::vector<std::complex<double>>>>
      Couplings_Lepton_Higgs_21;
 
  std::vector<std::vector<std::vector<double>>> LambdaGauge_3;
  std::vector<std::vector<std::vector<double>>> LambdaHiggs_3;
  std::vector<std::vector<std::vector<double>>> LambdaHiggs_3_CT;
  std::vector<std::vector<std::vector<std::complex<double>>>> LambdaQuark_3;
  std::vector<std::vector<std::vector<std::complex<double>>>> LambdaLepton_3;
  std::vector<std::vector<std::vector<std::vector<std::complex<double>>>>>
      LambdaQuark_4;
  std::vector<std::vector<std::vector<std::vector<std::complex<double>>>>>
      LambdaLepton_4;
 
  std::vector<std::vector<double>> DebyeHiggs;
  std::vector<std::vector<double>> DebyeGauge;
  std::vector<std::size_t> VevOrder;
 
public:
  [[deprecated("Will call Class_Potential_Origin with GetSMConstants(). "
               "Please use the "
               "detailed overload "
               "to ensure consistent SM constants through all "
               "routines.")]] Class_Potential_Origin();
 
  Class_Potential_Origin(const ISMConstants &smConstants);
  virtual ~Class_Potential_Origin();
 
  const auto &Get_Curvature_Higgs_L2() const { return Curvature_Higgs_L2; }
 
  const auto &Get_Curvature_Higgs_L3() const { return Curvature_Higgs_L3; }
 
  const auto &Get_Curvature_Higgs_L4() const { return Curvature_Higgs_L4; }
 
  const auto &Get_Curvature_Gauge_G2H2() const { return Curvature_Gauge_G2H2; }
 
  const auto &Get_Curvature_Lepton_F2H1() const
  {
    return Curvature_Lepton_F2H1;
  }
 
  const auto &Get_Curvature_Lepton_F2() const { return Curvature_Lepton_F2; }
 
  const auto &Get_Curvature_Quark_F2H1() const { return Curvature_Quark_F2H1; }
 
  const auto &Get_Curvature_Quark_F2() const { return Curvature_Quark_F2; }
 
  const std::vector<std::size_t> &Get_VevOrder() const { return VevOrder; }
 
  double get_scale() const { return scale; }
 
  void set_scale(double scale_new) { scale = scale_new; }
  std::size_t get_nPar() const { return nPar; }
  std::size_t get_nParCT() const { return nParCT; }
  std::size_t get_nVEV() const { return nVEV; }
  std::vector<double> get_vevTreeMin() const { return vevTreeMin; }
  double get_vevTreeMin(const std::size_t &k) const { return vevTreeMin.at(k); }
  std::vector<double> get_parStored() const { return parStored; }
  std::vector<double> get_parCTStored() const { return parCTStored; }
  std::size_t get_NGauge() const { return NGauge; }
  std::size_t get_NQuarks() const { return NQuarks; }
  std::size_t get_NHiggs() const { return NHiggs; }
  std::size_t get_NChargedHiggs() const { return NChargedHiggs; }
  std::size_t get_NNeutralHiggs() const { return NNeutralHiggs; }
  std::size_t get_NLepton() const { return NLepton; }
  ModelID::ModelIDs get_Model() const { return Model; }
 
  const std::vector<std::vector<double>> &get_DebyeHiggs() const
  {
    return DebyeHiggs;
  }
 
  void set_InputLineNumber(int InputLineNumber_in)
  {
    InputLineNumber = InputLineNumber_in;
  }
  double get_TripleHiggsCorrectionsTreePhysical(std::size_t i,
                                                std::size_t j,
                                                std::size_t k) const
  {
    return TripleHiggsCorrectionsTreePhysical.at(i).at(j).at(k);
  }
  double get_TripleHiggsCorrectionsCTPhysical(std::size_t i,
                                              std::size_t j,
                                              std::size_t k) const
  {
    return TripleHiggsCorrectionsCTPhysical.at(i).at(j).at(k);
  }
  double get_TripleHiggsCorrectionsCWPhysical(std::size_t i,
                                              std::size_t j,
                                              std::size_t k) const
  {
    return TripleHiggsCorrectionsCWPhysical.at(i).at(j).at(k);
  }
 
  double get_HiggsRotationMatrix(std::size_t i, std::size_t j) const
  {
    return HiggsRotationMatrix.at(i).at(j);
  }
  double get_HiggsRotationMatrixEnsuredConvention(std::size_t i,
                                                  std::size_t j) const
  {
    return HiggsRotationMatrixEnsuredConvention.at(i).at(j);
  }
 
  void setUseIndexCol(std::string legend);
 
  bool getUseIndexCol();
 
  void sym2Dim(std::vector<std::vector<double>> &Tensor2Dim,
               std::size_t Nk1,
               std::size_t Nk2);
 
  void sym3Dim(std::vector<std::vector<std::vector<double>>> &Tensor3Dim,
               std::size_t Nk1,
               std::size_t Nk2,
               std::size_t Nk3);
 
  void sym4Dim(
      std::vector<std::vector<std::vector<std::vector<double>>>> &Tensor4Dim,
      std::size_t Nk1,
      std::size_t Nk2,
      std::size_t Nk3,
      std::size_t Nk4);
 
  void initVectors();
 
  [[deprecated("Will call VEff with enum class Order. Consider using this "
               "instead.")]] virtual double
  VEff(const std::vector<double> &v, double Temp, int diff, int Order) const;
 
  virtual double VEff(const std::vector<double> &v,
                      double Temp        = 0,
                      int diff           = 0,
                      const Order &order = Order::OneLoop) const;
 
  double VTree(const std::vector<double> &v,
               int diff                      = 0,
               bool ForceExplicitCalculation = false) const;
  double CounterTerm(const std::vector<double> &v,
                     int diff                      = 0,
                     bool ForceExplicitCalculation = false) const;
  double V1Loop(const std::vector<double> &v, double Temp, int diff) const;
 
  virtual double EWSBVEV(const std::vector<double> &v) const;
 
  virtual void SetEWVEVZero(std::vector<double> &sol) const;
 
  virtual void ReadAndSet(const std::string &linestr,
                          std::vector<double> &par) = 0;
  virtual std::vector<std::string> addLegendCT() const = 0;
  virtual std::vector<std::string> addLegendTemp() const = 0;
  virtual std::vector<std::string> addLegendTripleCouplings() const = 0;
  virtual std::vector<std::string> addLegendVEV() const = 0;
 
  virtual void set_gen(const std::vector<double> &par) = 0;
  virtual void set_CT_Pot_Par(const std::vector<double> &par) = 0;
  virtual void write() const = 0;
 
  void set_All(const std::vector<double> &par,
               const std::vector<double> &parCT);
 
  virtual void SetCurvatureArrays() = 0;
  void CalculatePhysicalCouplings();
  std::vector<double> WeinbergFirstDerivative() const;
  std::vector<double> WeinbergSecondDerivative() const;
  Eigen::MatrixXd WeinbergSecondDerivativeAsMatrixXd() const;
  std::vector<double> WeinbergThirdDerivative() const;
  std::vector<double> WeinbergForthDerivative() const;
 
  void CalculateDebye(bool forceCalculation = false);
 
  void CalculateDebyeGauge();
  void Prepare_Triple();
 
  virtual bool CalculateDebyeSimplified() = 0;
 
  virtual bool CalculateDebyeGaugeSimplified() = 0;
 
  virtual double VTreeSimplified(const std::vector<double> &v) const = 0;
  bool UseVTreeSimplified = false;
  virtual double VCounterSimplified(const std::vector<double> &v) const = 0;
  bool UseVCounterSimplified = false;
 
  std::vector<double> HiggsMassesSquared(const std::vector<double> &v,
                                         const double &Temp = 0,
                                         const int &diff    = 0) const;
 
  Eigen::MatrixXd HiggsMassMatrix(const std::vector<double> &v,
                                  double Temp = 0,
                                  int diff    = 0) const;
 
  std::vector<double> GaugeMassesSquared(const std::vector<double> &v,
                                         const double &Temp = 0,
                                         const int &diff    = 0) const;
  std::vector<double> QuarkMassesSquared(const std::vector<double> &v,
                                         const int &diff = 0) const;
  std::vector<double> LeptonMassesSquared(const std::vector<double> &v,
                                          const int &diff = 0) const;
 
  std::vector<std::complex<double>>
  QuarkMasses(const std::vector<double> &v) const;
  Eigen::MatrixXcd QuarkMassMatrix(const std::vector<double> &v) const;
  std::vector<std::complex<double>>
  LeptonMasses(const std::vector<double> &v) const;
  Eigen::MatrixXcd LeptonMassMatrix(const std::vector<double> &v) const;
 
  const double ARMZeroThreshold = 1e-5;
  virtual void AdjustRotationMatrix() = 0;
  bool CheckRotationMatrix();
  virtual void TripleHiggsCouplings() = 0;
  double fbase(double MassSquaredA, double MassSquaredB) const;
  double
  fbaseTri(double MassSquaredA, double MassSquaredB, double MassSquaredC) const;
  double fbaseFour(double MassSquaredA,
                   double MassSquaredB,
                   double MassSquaredC,
                   double MassSquaredD) const;
 
  virtual std::vector<double> calc_CT() const = 0;
 
  double CWTerm(double MassSquared, double cb, int diff = 0) const;
  double FCW(double MassSquared) const;
  double boson(double MassSquared, double Temp, double cb, int diff = 0) const;
  [[deprecated(
      "Use this only if you want to calculate the effective potential with the "
      "exact same routine used in BSMPT v1.x. Use boson otherwise.")]] double
  boson_legacy(double MassSquared, double Temp, double cb, int diff = 0) const;
  double fermion(double MassSquared, double Temp, int diff = 0) const;
  [[deprecated(
      "Use this only if you want to calculate the effective potential with the "
      "exact same routine used in BSMPT v1.x. Use fermion otherwise.")]] double
  fermion_legacy(double Mass, double Temp, int diff = 0) const;
  [[deprecated("Use this only if you want to calculate the effective potential "
               "with the exact same routine used in BSMPT v1.x. "
               "Otherwise use ThermalFunctions::JInterpolatedHigh.")]] double
  Vl(double MassSquared, double Temp, int n, int diff = 0) const;
  [[deprecated(
      "Use this only if you want to calculate the effective potential with the "
      "exact same routine used in BSMPT v1.x. "
      "Otherwise use ThermalFunctions::JfermionInterpolatedLow.")]] double
  Vsf(double MassSquared, double Temp, int n, int diff = 0) const;
 
  [[deprecated(
      "Use this only if you want to calculate the effective potential with the "
      "exact same routine used in BSMPT v1.x. "
      "Otherwise use ThermalFunctions::JbosonInterpolatedLow.")]] double
  Vsb(double MassSquaredIn, double Temp, int n, int diff = 0) const;
 
  void resetbools();
 
  std::vector<double>
  MinimizeOrderVEV(const std::vector<double> &VEVminimizer) const;
 
  std::vector<double>
  FirstDerivativeOfEigenvalues(const Eigen::Ref<Eigen::MatrixXcd> M,
                               const Eigen::Ref<Eigen::MatrixXcd> MDiff) const;
  std::vector<double> SecondDerivativeOfEigenvaluesNonRepeated(
      const Eigen::Ref<Eigen::MatrixXd> M,
      const Eigen::Ref<Eigen::MatrixXd> MDiffX,
      const Eigen::Ref<Eigen::MatrixXd> MDiffY,
      const Eigen::Ref<Eigen::MatrixXd> MDiffXY) const;
 
  bool CheckNLOVEV(const std::vector<double> &v) const;
 
  virtual void Debugging(const std::vector<double> &input,
                         std::vector<double> &output) const = 0;
 
  std::vector<std::vector<double>> SignSymmetries;
  void FindSignSymmetries();
 
  void SetUseTreeLevel(bool val);
 
  std::pair<std::vector<double>, std::vector<double>>
  initModel(std::string linestr);
 
  std::vector<double> initModel(const std::vector<double> &par,
                                const bool &adjust_rot_matrix = true);
 
  std::vector<double> resetScale(const double &newScale);
 
  double CalculateRatioAlpha(const std::vector<double> &vev_symmetric,
                             const std::vector<double> &vev_broken,
                             const double &Temp) const;
 
  Eigen::VectorXd NablaVCT(const std::vector<double> &v) const;
  Eigen::MatrixXd HessianCT(const std::vector<double> &v) const;
 
  virtual std::vector<double> GetCTIdentities() const;
};
 
} // namespace BSMPT
#endif /* CLASSPOTENTIALORIGIN_H_ */