Skip to content

Commit

Permalink
fist implementation of b0 function with squared arguments
Browse files Browse the repository at this point in the history
  • Loading branch information
Alexander Voigt authored and Alexander Voigt committed Feb 28, 2018
1 parent 3a80cae commit 822be54
Show file tree
Hide file tree
Showing 3 changed files with 115 additions and 8 deletions.
50 changes: 49 additions & 1 deletion src/numerics.cpp
View file
Expand Up @@ -69,9 +69,10 @@ T divide_finite(T a, T b) noexcept {
double fB(const std::complex<double>& a) noexcept
{
using flexiblesusy::fast_log;
using std::abs;
const double x = a.real();

if (fabs(x) < EPSTOL)
if (abs(x) < EPSTOL)
return -1. - x + sqr(x) * 0.5;

if (is_close(x, 1., EPSTOL))
Expand Down Expand Up @@ -418,4 +419,51 @@ double a0(double m2, double q2) noexcept
return m2 * (1.0 - std::log(std::abs(m2 / q2)));
}

/**
* B0 function with squared arguments, from hep-ph/9606211.
* Note it returns the REAL PART ONLY.
*/
double b0(double p2, double m12, double m22, double q2) noexcept
{
using namespace softsusy;
using std::log;
using std::abs;

// protect against infrared divergence
if (is_close(p2, 0., EPSTOL) && is_close(m12, 0., EPSTOL)
&& is_close(m22, 0., EPSTOL))
return 0.;

double ans = 0.;
const double mMax2 = std::max(abs(m12), abs(m22));
const double pTest = abs(divide_finite(p2, mMax2));

if (pTest > 1e-10) {
const double s = p2 - m22 + m12;
const double s2 = sqr(s);
const std::complex<double> iEpsilon(0., EPSTOL * mMax2);
const std::complex<double> xPlus =
(s + sqrt(s2 - 4.*p2*(m12 - iEpsilon))) / (2.*p2);
const std::complex<double> xMinus =
(s - sqrt(s2 - 4.*p2*(m12 - iEpsilon))) / (2.*p2);

ans = -log(abs(p2 / q2)) - fB(xPlus) - fB(xMinus);
} else {
const double mMin2 = std::min(abs(m12), abs(m22));

if (is_close(m12, m22, EPSTOL)) {
ans = -log(abs(m12 / q2));
} else {
if (mMin2 < 1.e-30) {
ans = 1. - log(abs(mMax2 / q2));
} else {
ans = (m12*(1. - log(abs(m12/q2))) - m22*(1. - log(abs(m22/q2))))
/ (m12 - m22);
}
}
}

return ans;
}

} // namespace flexiblesusy
1 change: 1 addition & 0 deletions src/numerics.h
View file
Expand Up @@ -47,6 +47,7 @@ double b22bar(double p, double m1, double m2, double q) noexcept;
namespace flexiblesusy {

double a0(double m2, double q2) noexcept;
double b0(double p2, double m12, double m22, double q2) noexcept;

} // namespace flexiblesusy

Expand Down
72 changes: 65 additions & 7 deletions test/test_pv.cpp
View file
Expand Up @@ -107,7 +107,7 @@ BOOST_AUTO_TEST_CASE( test_ReA0 )
BOOST_CHECK_CLOSE_FRACTION(softsusy::a0(1e-5, scale), flexiblesusy::a0(1e-10, scale2), 1e-15);
}

BOOST_AUTO_TEST_CASE( test_ReB0 )
BOOST_AUTO_TEST_CASE( test_ReB0_limits )
{
BOOST_CHECK_EQUAL(ReB0(0., 0., 0., scale2), 0.);

Expand All @@ -120,6 +120,70 @@ BOOST_AUTO_TEST_CASE( test_ReB0 )
BOOST_CHECK_EQUAL(ReB0(p2, 0., p2, scale2), ReB0(p2, p2, 0., scale2));
}

struct Values {
Values(double p_, double m1_, double m2_, double q_)
: p(p_), m1(m1_), m2(m2_), q(q_) {}
double p{}, m1{}, m2{}, q{};
};

BOOST_AUTO_TEST_CASE( test_ReB0_values )
{
const std::vector<Values> vals = {
Values(0. , 1. , 0., 1.),
Values(1e-1, 1. , 0., 1.),
Values(1e-2, 1. , 0., 1.),
Values(1e-3, 1. , 0., 1.),
Values(1e-4, 1. , 0., 1.),
Values(1e-5, 1. , 0., 1.),
Values(0. , 1. , 1e-15, 1.),
Values(1e-1, 1. , 1e-15, 1.),
Values(1e-2, 1. , 1e-15, 1.),
Values(1e-3, 1. , 1e-15, 1.),
Values(1e-4, 1. , 1e-15, 1.),
Values(1e-5, 1. , 1e-15, 1.),
Values(0. , 1e20, 0., 1.),
Values(1e-1, 1e20, 0., 1.),
Values(1e-2, 1e20, 0., 1.),
Values(1e-3, 1e20, 0., 1.),
Values(1e-4, 1e20, 0., 1.),
Values(1e-5, 1e20, 0., 1.),
Values(0. , 0. , 1., 1.),
Values(1e-1, 0. , 1., 1.),
Values(1e-2, 0. , 1., 1.),
Values(1e-3, 0. , 1., 1.),
Values(1e-4, 0. , 1., 1.),
Values(1e-5, 0. , 1., 1.),
Values(0. , 1e20, 1., 1.),
Values(1e-1, 1e20, 1., 1.),
Values(1e-2, 1e20, 1., 1.),
Values(1e-3, 1e20, 1., 1.),
Values(1e-4, 1e20, 1., 1.),
Values(1e-5, 1e20, 1., 1.),
Values(0. , 0., 1e20, 1.),
Values(1e-1, 0., 1e20, 1.),
Values(1e-2, 0., 1e20, 1.),
Values(1e-3, 0., 1e20, 1.),
Values(1e-4, 0., 1e20, 1.),
Values(1e-5, 0., 1e20, 1.),
Values(1. , 1. , 1., 1.),
Values(1. , 2. , 3., 4.)
};

for (const auto v: vals) {
const auto p = v.p;
const auto m1 = v.m1;
const auto m2 = v.m2;
const auto q = v.q;

const auto v1 = softsusy::b0(p,m1,m2,q);
const auto v2 = flexiblesusy::b0(sqr(p),sqr(m1),sqr(m2),sqr(q));

BOOST_TEST_MESSAGE("testing p = " << p << ", m1 = " << m1
<< ", m2 = " << m2 << ", q = " << q);
BOOST_CHECK_CLOSE_FRACTION(v1, v2, 8e-5);
}
}

BOOST_AUTO_TEST_CASE( test_ReB1 )
{
BOOST_CHECK_EQUAL(ReB1(0., 0., 0., scale2), 0.);
Expand All @@ -130,12 +194,6 @@ BOOST_AUTO_TEST_CASE( test_ReB1 )
BOOST_CHECK(ReB1(0,1e300,0,1e-10) != 0.);
}

struct Values {
Values(double p_, double m1_, double m2_, double q_)
: p(p_), m1(m1_), m2(m2_), q(q_) {}
double p{}, m1{}, m2{}, q{};
};

BOOST_AUTO_TEST_CASE( test_ReB1_integral )
{
const std::vector<Values> vals = {
Expand Down

0 comments on commit 822be54

Please sign in to comment.