Skip to content

Commit

Permalink
use new loop functions with squared arguments
Browse files Browse the repository at this point in the history
  • Loading branch information
Alexander Voigt authored and Alexander Voigt committed Mar 3, 2018
1 parent cdf6386 commit 4930d4b
Show file tree
Hide file tree
Showing 5 changed files with 76 additions and 44 deletions.
10 changes: 5 additions & 5 deletions src/pv.cpp
View file
Expand Up @@ -27,7 +27,7 @@
#elif defined(ENABLE_FFLITE)
# include "fflite.hpp"
#else
# include "numerics.h"
# include "pv2.hpp"
#endif

namespace flexiblesusy {
Expand Down Expand Up @@ -442,7 +442,7 @@ double ReA0(double m2, double scl2) noexcept
#if defined(ENABLE_LOOPTOOLS) || defined(ENABLE_FFLITE)
return A0(m2, scl2).real();
#else
return softsusy::a0(sqrt(m2), sqrt(scl2));
return a0(m2, scl2);
#endif
}

Expand All @@ -451,7 +451,7 @@ double ReB0(double p2, double m2a, double m2b, double scl2) noexcept
#if defined(ENABLE_LOOPTOOLS) || defined(ENABLE_FFLITE)
return B0(p2, m2a, m2b, scl2).real();
#else
return softsusy::b0(sqrt(p2), sqrt(m2a), sqrt(m2b), sqrt(scl2));
return b0(p2, m2a, m2b, scl2);
#endif
}

Expand All @@ -460,7 +460,7 @@ double ReB1(double p2, double m2a, double m2b, double scl2) noexcept
#if defined(ENABLE_LOOPTOOLS) || defined(ENABLE_FFLITE)
return B1(p2, m2a, m2b, scl2).real();
#else
return -softsusy::b1(sqrt(p2), sqrt(m2a), sqrt(m2b), sqrt(scl2));
return -b1(p2, m2a, m2b, scl2);
#endif
}

Expand All @@ -469,7 +469,7 @@ double ReB00(double p2, double m2a, double m2b, double scl2) noexcept
#if defined(ENABLE_LOOPTOOLS) || defined(ENABLE_FFLITE)
return B00(p2, m2a, m2b, scl2).real();
#else
return softsusy::b22(sqrt(p2), sqrt(m2a), sqrt(m2b), sqrt(scl2));
return b22(p2, m2a, m2b, scl2);
#endif
}

Expand Down
2 changes: 1 addition & 1 deletion src/pv.hpp
View file
Expand Up @@ -99,7 +99,7 @@ std::complex<double> H0(T p2, T m2a, T m2b, double scl2) noexcept
#endif

// the following are mainly for interfacing with loop function
// implementations from softsusy since they come only with double
// implementations inspired by softsusy since they come only with double
// return type. If LoopTools or FF is in use, they reduce simply to
// A0(m2, scl2).real(), etc.
double ReA0 (double m2, double scl2) PVATTR;
Expand Down
6 changes: 3 additions & 3 deletions src/threshold_loop_functions.cpp
View file
Expand Up @@ -20,7 +20,7 @@
#include "dilog.hpp"
#include "pv.hpp"
#include "logger.hpp"
#include "numerics.h"
#include "pv2.hpp"

#include <cmath>
#include <limits>
Expand Down Expand Up @@ -1822,12 +1822,12 @@ double DB0(double m1, double m2)
*/
double C0(double m1, double m2, double m3)
{
return softsusy::c0(m1, m2, m3);
return c0(m1*m1, m2*m2, m3*m3);
}

double D0(double m1, double m2, double m3, double m4)
{
return softsusy::d0(m1,m2,m3,m4);
return d0(m1*m1,m2*m2,m3*m3,m4*m4);
}

/**
Expand Down
88 changes: 60 additions & 28 deletions src/weinberg_angle.cpp
View file
Expand Up @@ -22,7 +22,7 @@
#include "logger.hpp"
#include "numerics2.hpp"
#include "config.h"
#include "numerics.h"
#include "pv2.hpp"

#include <limits>
#include <cmath>
Expand All @@ -41,8 +41,6 @@ namespace flexiblesusy {

namespace weinberg_angle {

using namespace softsusy;

Weinberg_angle::Data::Data()
: scale(0.)
, alpha_em_drbar(0.)
Expand Down Expand Up @@ -236,6 +234,40 @@ int Weinberg_angle::calculate(double rho_start, double sin_start)
return no_convergence_error;
}

namespace {

double B0(double p, double m1, double m2, double q) noexcept
{
return b0(p*p, m1*m1, m2*m2, q*q);
}

double B1(double p, double m1, double m2, double q) noexcept
{
return b1(p*p, m1*m1, m2*m2, q*q);
}

double H0(double p, double m1, double m2, double q) noexcept
{
return h0(p*p, m1*m1, m2*m2, q*q);
}

double C0(double m1, double m2, double m3) noexcept
{
return c0(m1*m1, m2*m2, m3*m3);
}

double D0(double m1, double m2, double m3, double m4) noexcept
{
return d0(m1*m1, m2*m2, m3*m3, m4*m4);
}

double D27(double m1, double m2, double m3, double m4) noexcept
{
return d27(m1*m1, m2*m2, m3*m3, m4*m4);
}

} // anonymous namespace

/**
* Calculates the \f$\Delta\hat{\rho}\f$ corrections as defined in
* Eqs. (C.4), (C.6) from hep-ph/9606211 .
Expand Down Expand Up @@ -546,12 +578,12 @@ double Weinberg_angle::calculate_delta_vb_susy(

double deltaZnue = 0.0, deltaZe = 0.0;
for (int i = 0; i < dimN; i++) {
deltaZnue += - Sqr(Abs(bChi0NuNul(i))) * b1(0.0, mneut(i), msnue, q);
deltaZe += - Sqr(Abs(bChi0ESell(i))) * b1(0.0, mneut(i), mselL, q);
deltaZnue += - Sqr(Abs(bChi0NuNul(i))) * B1(0.0, mneut(i), msnue, q);
deltaZe += - Sqr(Abs(bChi0ESell(i))) * B1(0.0, mneut(i), mselL, q);
}
for (int i = 0; i < dimC; i++) {
deltaZnue += - Sqr(Abs(bChicNuSell(i))) * b1(0.0, mch(i), mselL, q);
deltaZe += - Sqr(Abs(aChicESnul(i))) * b1(0.0, mch(i), msnue, q);
deltaZnue += - Sqr(Abs(bChicNuSell(i))) * B1(0.0, mch(i), mselL, q);
deltaZe += - Sqr(Abs(aChicESnul(i))) * B1(0.0, mch(i), msnue, q);
}

Eigen::VectorXd bPsicNuSmul(Eigen::VectorXd::Zero(dimC));
Expand All @@ -571,12 +603,12 @@ double Weinberg_angle::calculate_delta_vb_susy(

double deltaZnumu = 0.0, deltaZmu = 0.0;
for (int i = 0; i < dimN; i++) {
deltaZnumu += - Sqr(Abs(bChi0NuNul(i))) * b1(0.0, mneut(i), msnumu, q);
deltaZmu += - Sqr(Abs(bChi0MuSmul(i))) * b1(0.0, mneut(i), msmuL, q);
deltaZnumu += - Sqr(Abs(bChi0NuNul(i))) * B1(0.0, mneut(i), msnumu, q);
deltaZmu += - Sqr(Abs(bChi0MuSmul(i))) * B1(0.0, mneut(i), msmuL, q);
}
for (int i = 0; i < dimC; i++) {
deltaZnumu += - Sqr(Abs(bChicNuSmul(i))) * b1(0.0, mch(i), msmuL, q);
deltaZmu += - Sqr(Abs(aChicMuSnul(i))) * b1(0.0, mch(i), msnumu, q);
deltaZnumu += - Sqr(Abs(bChicNuSmul(i))) * B1(0.0, mch(i), msmuL, q);
deltaZmu += - Sqr(Abs(aChicMuSnul(i))) * B1(0.0, mch(i), msnumu, q);
}

Eigen::MatrixXd aPsi0PsicW(Eigen::MatrixXd::Zero(dimN,dimC)),
Expand All @@ -593,15 +625,15 @@ double Weinberg_angle::calculate_delta_vb_susy(
aChi0ChicW = n.conjugate() * aPsi0PsicW * v.transpose();
bChi0ChicW = n * bPsi0PsicW * u.adjoint();

const double b0_0_mselL_msnue = b0(0.0, mselL, msnue, q);
const double b0_0_mnmuL_msnumu = b0(0.0, msmuL, msnumu, q);
const double b0_0_mselL_msnue = B0(0.0, mselL, msnue, q);
const double b0_0_mnmuL_msnumu = B0(0.0, msmuL, msnumu, q);

std::complex<double> deltaVE;
for (int i = 0; i < dimC; i++) {
for (int j = 0; j < dimN; j++) {
const double c0_mselL_mch_mneut = c0(mselL, mch(i), mneut(j));
const double c0_msnue_mch_mneut = c0(msnue, mch(i), mneut(j));
const double b0_0_mch_mneut = b0(0.0, mch(i), mneut(j), q);
const double c0_mselL_mch_mneut = C0(mselL, mch(i), mneut(j));
const double c0_msnue_mch_mneut = C0(msnue, mch(i), mneut(j));
const double b0_0_mch_mneut = B0(0.0, mch(i), mneut(j), q);

deltaVE += bChicNuSell(i) * Conj(bChi0ESell(j)) *
(- ROOT2 / g * aChi0ChicW(j, i) * mch(i) * mneut(j) *
Expand All @@ -622,15 +654,15 @@ double Weinberg_angle::calculate_delta_vb_susy(
deltaVE +=
0.5 * Conj(bChi0ESell(j)) * bChi0NuNul(j) *
(b0_0_mselL_msnue + Sqr(mneut(j)) *
c0(mneut(j), mselL, msnue) + 0.5);
C0(mneut(j), mselL, msnue) + 0.5);
}

std::complex<double> deltaVMu;
for (int i = 0; i < dimC; i++) {
for (int j = 0; j < dimN; j++) {
const double c0_msmuL_mch_mneut = c0(msmuL, mch(i), mneut(j));
const double c0_msnumu_mch_mneut = c0(msnumu, mch(i), mneut(j));
const double b0_0_mch_mneut = b0(0.0, mch(i), mneut(j), q);
const double c0_msmuL_mch_mneut = C0(msmuL, mch(i), mneut(j));
const double c0_msnumu_mch_mneut = C0(msnumu, mch(i), mneut(j));
const double b0_0_mch_mneut = B0(0.0, mch(i), mneut(j), q);

deltaVMu += bChicNuSmul(i) * Conj(bChi0MuSmul(j)) *
(- ROOT2 / g * aChi0ChicW(j, i) * mch(i) * mneut(j) *
Expand All @@ -651,24 +683,24 @@ double Weinberg_angle::calculate_delta_vb_susy(
deltaVMu +=
0.5 * Conj(bChi0MuSmul(j)) * bChi0NuNul(j) *
(b0_0_mnmuL_msnumu + Sqr(mneut(j)) *
c0(mneut(j), msmuL, msnumu) + 0.5);
C0(mneut(j), msmuL, msnumu) + 0.5);
}

std::complex<double> a1;
for(int i = 0; i < dimC; i++) {
for(int j = 0; j < dimN; j++) {
a1 += 0.5 * aChicMuSnul(i) * Conj(bChicNuSell(i)) *
bChi0NuNul(j) * bChi0ESell(j) * mch(i) * mneut(j) *
d0(mselL, msnumu, mch(i), mneut(j));
D0(mselL, msnumu, mch(i), mneut(j));
a1 += 0.5 * Conj(aChicESnul(i)) * bChicNuSmul(i) *
Conj(bChi0NuNul(j)) * Conj(bChi0MuSmul(j)) * mch(i) * mneut(j) *
d0(msmuL, msnue, mch(i), mneut(j));
D0(msmuL, msnue, mch(i), mneut(j));
a1 += bChicNuSmul(i) * Conj(bChicNuSell(i)) *
Conj(bChi0MuSmul(j)) * bChi0ESell(j) *
d27(msmuL, mselL, mch(i), mneut(j));
D27(msmuL, mselL, mch(i), mneut(j));
a1 += Conj(aChicMuSnul(i)) * aChicESnul(i) *
bChi0NuNul(j) * Conj(bChi0NuNul(j)) *
d27(msnumu, msnue, mch(i), mneut(j));
D27(msnumu, msnue, mch(i), mneut(j));
}
}

Expand Down Expand Up @@ -738,7 +770,7 @@ double Weinberg_angle::calculate_self_energy_w_top(
const double g2 = data.g2;

const double self_energy_w =
0.5 * Nc * hfn(p, mt, mb, q) * Sqr(g2) * oneOver16PiSqr;
0.5 * Nc * H0(p, mt, mb, q) * Sqr(g2) * oneOver16PiSqr;

return self_energy_w;
}
Expand Down Expand Up @@ -767,8 +799,8 @@ double Weinberg_angle::calculate_self_energy_z_top(
const double guR = 2.0 * sw2 / 3.0;

const double self_energy_z =
Nc * (hfn(p, mt, mt, q) * (Sqr(guL) + Sqr(guR))
- 4.0 * guL * guR * Sqr(mt) * b0(p, mt, mt, q))
Nc * (H0(p, mt, mt, q) * (Sqr(guL) + Sqr(guR))
- 4.0 * guL * guR * Sqr(mt) * B0(p, mt, mt, q))
* Sqr(g2) / cw2 * oneOver16PiSqr;

return self_energy_z;
Expand Down
14 changes: 7 additions & 7 deletions templates/weinberg_angle.cpp.in
View file
Expand Up @@ -25,9 +25,9 @@
#include "logger.hpp"
#include "numerics2.hpp"
#include "config.h"
#include "numerics.h"
#include "error.hpp"
#include "pv.hpp"
#include "pv2.hpp"

#include <limits>
#include <cmath>
Expand Down Expand Up @@ -420,17 +420,17 @@ double CLASSNAME::B1(double p2, double m12, double m22) const noexcept

double CLASSNAME::C0(double m12, double m22, double m32) const noexcept
{
return softsusy::c0(std::sqrt(m12), std::sqrt(m22), std::sqrt(m32));
return c0(m12, m22, m32);
}

double CLASSNAME::D0(double m12, double m22, double m32, double m42) const noexcept
{
return softsusy::d0(std::sqrt(m12), std::sqrt(m22), std::sqrt(m32), std::sqrt(m42));
return d0(m12, m22, m32, m42);
}

double CLASSNAME::D27(double m12, double m22, double m32, double m42) const noexcept
{
return softsusy::d27(std::sqrt(m12), std::sqrt(m22), std::sqrt(m32), std::sqrt(m42));
return d27(m12, m22, m32, m42);
}

/**
Expand Down Expand Up @@ -546,8 +546,8 @@ double CLASSNAME::calculate_self_energy_@VectorZ@_top(double p, double mt) const

const double self_energy_z_top =
Nc * Sqr(g2) / cw2 * oneOver16PiSqr *
(softsusy::hfn(p, mt, mt, q) * (Sqr(guL) + Sqr(guR)) -
4.0 * guL * guR * Sqr(mt) * softsusy::b0(p, mt, mt, q));
(h0(p*p, mt*mt, mt*mt, q*q) * (Sqr(guL) + Sqr(guR)) -
4.0 * guL * guR * Sqr(mt) * b0(p*p, mt*mt, mt*mt, q*q));

return self_energy_z_top;
}
Expand All @@ -568,7 +568,7 @@ double CLASSNAME::calculate_self_energy_@VectorW@_top(double p, double mt) const
@DefSMleftCoupling@

const double self_energy_w_top =
0.5 * Nc * softsusy::hfn(p, mt, mb, q) * Sqr(g2) * oneOver16PiSqr;
0.5 * Nc * h0(p*p, mt*mt, mb*mb, q*q) * Sqr(g2) * oneOver16PiSqr;

return self_energy_w_top;
}
Expand Down

0 comments on commit 4930d4b

Please sign in to comment.