$ | \Delta B = 2 | $ Amplitude Class More...

#include <AmpDB2.h>

Detailed Description

$ | \Delta B = 2 | $ Amplitude Class

Author: HEPfit Collaboration

Copyright: GNU General Public License

This class is related to the calculation of the $ B_{d,s}-\bar{B}_{d,s}$ mixing.

Definition at line 36 of file AmpDB2.h.

Public Types
enum	mass_schemes { pole , MSbar , PS , only1overmb , MSbar_partialNNLO , PS_partialNNLO , MSbar_partialN3LO , PS_partialN3LO }

enum	quark { d , s }

enum	quarks { cc , cu , uu }

Public Member Functions
	AmpDB2 (const StandardModel &SM_i, int BMeson_i, bool flag_fixmub=false, bool flag_RI=false)
	Constructor. More...

gslpp::complex	getGamma21overM21 (orders order, mass_schemes mass_scheme=MSbar)
	The value of $\frac{\Gamma_{21},M_{21}}$ from Gerlach (2205.07907 and thesis) More...

gslpp::complex	getGamma21overM21_tradBasis (orders order)
	The value of $\frac{\Gamma_{21},M_{21}}$ in the traditional basis. More...

gslpp::complex	getM21 (orders order)
	The value of $M_{21}$ for $B_{d,s}$ mesons. More...

gslpp::complex	getPB ()

gslpp::complex	getRB (orders order)

Protected Member Functions
gslpp::complex	Gamma21overM21 (orders order, mass_schemes mass_scheme, int BMeson)
	A method to compute $\frac{\Gamma_{21},M_{21}}^{bq}$ @detail source: Marvin Gerlach (2205.07907 and thesis) with 1/mb corrections from Lenz (hep-ph/0612167) More...

gslpp::complex	Gamma21overM21_tradBasis (orders order, quark q)
	A method to compute $\frac{\Gamma_{21},M_{21}}^{bq}$ in the traditional basis @detail source: Ciuchini (hep-ph/0308029v2) More...

gslpp::complex	M21_Bd (orders order)
	A method to compute $M_{21}^{bd}$. More...

gslpp::complex	M21_Bs (orders order)
	A method to compute $M_{21}^{bs}$. More...

gslpp::complex	PBd ()

gslpp::complex	PBs ()

gslpp::complex	RBd (orders order)
	A method to compute the ratio of the absolute value of the $B_d$ mixing amplitude over the Standard Model value. More...

gslpp::complex	RBs (orders order)
	A method to compute the ratio of the absolute value of the $B_s$ mixing amplitude over the Standard Model value. More...

Private Member Functions
gslpp::vector< gslpp::complex >	c (quark q, orders order)
	Values of DB=2 Wilson coefficients from (hep-ph/0308029v2) transformed to the new basis. More...

gslpp::vector< gslpp::complex >	c_H (quark q, orders order)

gslpp::vector< gslpp::complex >	c_H_partial (quark q, int i)

void	compute_deltas_1overm (quark q)

void	compute_deltas_1overm_tradBasis (quark q)

void	compute_g ()

void	compute_matrixelements (quark q, orders order)
	A method to compute all DB=2 Wilson coefficients (me, me_R, me_Rtilde) More...

void	compute_partialNNLO ()

void	compute_pp_s ()

void	computeCKMandMasses (orders order=NNLO, mass_schemes mass_scheme=MSbar)
	A method to compute CKM elements, quark masses and alpha_s. More...

void	computeD (orders order)

void	computeF0 ()

void	computeF1 ()

void	computeP ()

void	computeWilsonCoeffsBuras ()

void	computeWilsonCoeffsMisiak ()

gslpp::complex	D (quarks qq, int k)

gslpp::complex	delta_1overm (quark q)
	Value of 1/mb corrections of $\Gamma_{21}$ (hep-ph/0612167) More...

gslpp::complex	delta_1overm_tradBasis (quark q)
	Value of 1/mb corrections of $\Gamma_{21}$ (hep-ph/0308029v2) More...

gslpp::complex	deltas_1overm (quarks qq, quark q)

gslpp::complex	deltas_1overm_tradBasis (quarks qq, quark q)

double	F (quarks qq, int k, int i, int j)

double	F0 (quarks qq, int k, int i, int j)

double	F1 (quarks qq, int k, int i, int j)

gslpp::complex	g (quarks qq, int i)

gslpp::complex	gtilde (quarks qq, int i)

gslpp::complex	H (quarks qq, orders order)

gslpp::vector< gslpp::complex >	H_allpartial (quarks qq)

gslpp::complex	H_partial (quarks qq, int i_start, int i_end, int j_start, int j_end, int n)

gslpp::complex	H_s (quarks qq, orders order)

gslpp::vector< gslpp::complex >	H_s_allpartial (quarks qq)

gslpp::complex	H_s_partial (quarks qq, int i_start, int i_end, int j_start, int j_end, int n)

int	index_deltas (quarks qq, quark q)

int	index_p (quarks qq, int i, int j, int n)

int	indexD (quarks qq, int k)

int	indexF (quarks qq, int k, int i, int j)

int	indexg (quarks qq, int i)

int	indexP (quarks qq, int k, int i, int j)

double	p (quarks qq, int i, int j, int n, bool flag_LOz=false)

double	P (quarks qq, int k, int i, int j)

double	p_s (quarks qq, int i, int j, int n, bool flag_LOz=false)

void	poletoMSbar_pp_s ()

void	poletoMSbar_pp_s_partialN3LO ()

void	poletoPS_pp_s ()

void	poletoPS_pp_s_partialN3LO ()

gslpp::vector< gslpp::complex >	transformation (gslpp::vector< gslpp::complex > result, orders order)

Private Attributes
double	a1 = 31./3. - 10./9. * nl

double	a2

double	as_4pi

double	as_4pi_mu1

double	as_4pi_mu2

double	b0 = 11. - 2. * nl/3.

double	b0_2 = b0 * b0

double	b1 = 102. - 38. * nl/3.

int	BMeson

gslpp::complex	C1_LO_1overm

gslpp::complex	C2_LO_1overm

gslpp::complex	C_1_SM

gslpp::complex	C_Buras_LO [8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

gslpp::complex	C_Buras_NLO [8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

gslpp::complex	C_Buras_NNLO [8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

gslpp::complex	C_Misiak_LO [8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

gslpp::complex	C_Misiak_NLO [8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

gslpp::complex	C_Misiak_NNLO [8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

gslpp::complex	cache_deltas_1overm [6] = {0.}

gslpp::complex	cache_deltas_1overm_tradBasis [6] = {0.}

double	cache_p [768] = { 0. }

double	cache_p_LO [576] = {0.}

double	cache_ps [768] = { 0. }

double	cache_ps_LO [576] = {0.}

gslpp::complex	cacheC [8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

gslpp::complex	cacheD [6] = {0.}

double	cacheF0 [24] = {0.}

double	cacheF1 [24] = {0.}

gslpp::complex	cacheg [12] = {0.}

gslpp::complex	cachegtilde [12] = {0.}

double	cacheP [84] = {0.}

const double	Cl2PI3 = 1.014941606

gslpp::matrix< double >	coeffsMStoRI

double	Dilogsigma

double	Dilogsigma2

double	Dilogz

bool	flag_fixmub

bool	flag_LOz = true

bool	flag_resumz

bool	flag_RI

double	Gf2

gslpp::complex	K_1

gslpp::complex	K_2

gslpp::complex	lambda_c_d

gslpp::complex	lambda_c_s

gslpp::complex	lambda_u_d

gslpp::complex	lambda_u_s

const double	log12sqrt52 = log(0.5 + sqrt5 / 2.)

double	log1minus4z

double	log1minusz

const double	log2 = log(2)

const double	log2_2 = log2 * log2

const double	log2_4 = log2_2 * log2_2

double	log2sigma

double	log2z

const double	log3 = log(3)

double	logsigma

double	logx_1

double	logx_2

double	logz

const double	M_PI2 = M_PI * M_PI

const double	M_PI3 = M_PI2 * M_PI

const double	M_PI4 = M_PI2 * M_PI2

double	Mb

double	MB

double	Mb2_prefactor

double	Mb2_prefactor_1overm

double	Mb_Mb

double	Mb_pole

double	Mb_PS

double	MB_s

double	Mc

double	Md

gslpp::vector< double >	me = gslpp::vector<double>(5, 0.)

gslpp::vector< double >	me_R = gslpp::vector<double>(5, 0.)

gslpp::vector< double >	me_Rtilde = gslpp::vector<double>(3, 0.)

gslpp::matrix< double >	meMStoRI

gslpp::vector< double >	meoverme0 = gslpp::vector<double>(3, 0.)

double	Ms

double	mu_1

double	mu_1_overm

double	mu_2

double	mu_b

double	mu_f = 2.

const StandardModel &	mySM

double	nl = 4.

double	oneminusz2

double	oneminusz_1overm2

double	PoletoMS_as1

double	PoletoMS_as2

double	PoletoMS_as2_z0

double	PoletoMS_as2_z1

double	PoletoMS_as3

double	PoletoPS_as1

double	PoletoPS_as2

double	PoletoPS_as3

const double	polylog4_12 = 0.517479062

double	sigma

double	sqrt1minus4z

double	sqrt1minus4z_1overm

const double	sqrt3 = sqrt(3)

const double	sqrt5 = sqrt(5)

double	sqrtz

const double	t_2 = -0.389011713

gslpp::complex	VcbVcd

gslpp::complex	VcbVcd2

gslpp::complex	VcbVcs

gslpp::complex	VcbVcs2

gslpp::complex	VtbVtd

gslpp::complex	VtbVtd2

gslpp::complex	VtbVts

gslpp::complex	VtbVts2

double	x_1

double	x_2

double	z

double	z2

double	z3

double	z4

double	z_1overm

double	z_1overm2

const double	zeta2 = gslpp_special_functions::zeta(2)

const double	zeta3 = gslpp_special_functions::zeta(3)

const double	zeta4 = gslpp_special_functions::zeta(4)

const double	zeta5 = gslpp_special_functions::zeta(5)

Member Enumeration Documentation

◆ mass_schemes

enum AmpDB2::mass_schemes

Enumerator
pole
MSbar
PS
only1overmb
MSbar_partialNNLO
PS_partialNNLO
MSbar_partialN3LO
PS_partialN3LO

Definition at line 63 of file AmpDB2.h.

63{pole, MSbar, PS, only1overmb, MSbar_partialNNLO, PS_partialNNLO, MSbar_partialN3LO, PS_partialN3LO}; //mass schemes

AmpDB2::MSbar_partialN3LO

@ MSbar_partialN3LO

Definition: AmpDB2.h:63

AmpDB2::MSbar

@ MSbar

Definition: AmpDB2.h:63

AmpDB2::PS

@ PS

Definition: AmpDB2.h:63

AmpDB2::PS_partialN3LO

@ PS_partialN3LO

Definition: AmpDB2.h:63

AmpDB2::pole

@ pole

Definition: AmpDB2.h:63

AmpDB2::MSbar_partialNNLO

@ MSbar_partialNNLO

Definition: AmpDB2.h:63

AmpDB2::only1overmb

@ only1overmb

Definition: AmpDB2.h:63

AmpDB2::PS_partialNNLO

@ PS_partialNNLO

Definition: AmpDB2.h:63

◆ quark

enum AmpDB2::quark

Enumerator
d
s

Definition at line 64 of file AmpDB2.h.

64{d,s}; /*quark index i used for $B_i$*/

AmpDB2::d

@ d

Definition: AmpDB2.h:64

s

Test Observable.

◆ quarks

enum AmpDB2::quarks

Enumerator
cc
cu
uu

Definition at line 65 of file AmpDB2.h.

65{cc, cu, uu}; /*combinations of u- and c- quarks in diagrams */

AmpDB2::cc

@ cc

Definition: AmpDB2.h:65

AmpDB2::cu

@ cu

Definition: AmpDB2.h:65

AmpDB2::uu

@ uu

Definition: AmpDB2.h:65

Constructor & Destructor Documentation

◆ AmpDB2()

AmpDB2::AmpDB2	(	const StandardModel &	SM_i,
		int	BMeson_i,
		bool	flag_fixmub = `false`,
		bool	flag_RI = `false`
	)

Constructor.

Parameters

[in]	SM_i	a reference to an object of type StandardModel
[in]	BMeson_i	an integer to specify the B meson: 0 for Bd and 1 for Bs

Definition at line 13 of file AmpDB2.cpp.

    : mySM(SM_i), meMStoRI(5, 0.), coeffsMStoRI(3, 0.)
{
    BMeson = BMeson_i;
    this->flag_resumz = true;
    this->flag_fixmub = flag_fixmub;
    this->flag_RI = flag_RI;
    switch (BMeson)
    {
    case 0:
        mySM.initializeBParameter("BBd");
        mySM.initializeBParameter("BBd_subleading");   
        mySM.initializeMeson(QCD::B_D);
        break;
    case 1:
        mySM.initializeBParameter("BBs");
        mySM.initializeBParameter("BBs_subleading");
        mySM.initializeMeson(QCD::B_S);
        break;
    default:
        throw std::runtime_error("AmpDB2::AmpDB2(): invalid B meson index");
        break;
    }
 
    // hep-ph/0606197 eq. 4.7 - 4.10
    double meMStoRI0[5] = {-3. - 5. / 3. + 8. * log2, 0., 0., 0., 0.},
           meMStoRI1[5] = {0., 61. / 9. + 44. / 9. * log2, -7. / 9. + 28. / 9. * log2, 0., 0.},
           meMStoRI2[5] = {0., -25. / 9. + 28. / 9. * log2, -29. / 9. + 44. / 9. * log2, 0., 0.},
           meMStoRI3[5] = {0., 0., 0., -5. / 3. + 13. - 2. / 3. * log2, -3. + 1. + 2. * log2},
           meMStoRI4[5] = {0., 0., 0., -7. / 2. + 11. / 2. + 2. * log2, -1. / 6. - 1. / 2. - 2. / 3. * log2};
    meMStoRI.assign(0, meMStoRI0);
    meMStoRI.assign(1, meMStoRI1);
    meMStoRI.assign(2, meMStoRI2);
    meMStoRI.assign(3, meMStoRI3);
    meMStoRI.assign(4, meMStoRI4);
    for (int i = 0; i <= 2; i++)
    {
        for (int j = 0; j <= 2; j++)
        {
            coeffsMStoRI.assign(i, j, meMStoRI(i, j));
        }
    }
}

Member Function Documentation

◆ c()

gslpp::vector< gslpp::complex > AmpDB2::c	(	quark	q,
		orders	order
	)

private

Values of DB=2 Wilson coefficients from (hep-ph/0308029v2) transformed to the new basis.

Parameters

[in]	quark	index of the neutral B mesons
[in]	order	the QCD order for the transformations @detail requires computeCKMandMasses() and "D"

Definition at line 875 of file AmpDB2.cpp.

{
    gslpp::vector<gslpp::complex> result = gslpp::vector<gslpp::complex>(3, 0.);
    switch (q)
    {
    case d:
        for (int i = 1; i <= 2; i++)
        {
            result.assign(i - 1,
                          VtbVtd2 * D(uu, i) + 2. * VcbVcd * VtbVtd * (D(uu, i) - D(cu, i)) + VcbVcd2 * (D(cc, i) + D(uu, i) - 2. * D(cu, i)));
        }
        break;
    case s:
        for (int i = 1; i <= 2; i++)
        {
            result.assign(i - 1,
                          VtbVts2 * D(uu, i) + 2. * VcbVcs * VtbVts * (D(uu, i) - D(cu, i)) + VcbVcs2 * (D(cc, i) + D(uu, i) - 2. * D(cu, i)));
        }
        break;
    default:
        throw std::runtime_error("AmpDB2::c(quark q, double mu_2): invalid quark index: ");
    }
 
    // transformation to the new basis (hep-ph/0612167 eq. 18 + 21)
    if (order == LO)
    {
        result.assign(0, result(0) - 0.5 * result(1));
        result.assign(2, -result(1));
        result.assign(1, 0.);
    }
    else if (order == NLO)
    {
        // save NLO contribution from transformation of the tradBasis to RI scheme
        gslpp::vector<gslpp::complex> RI_NLO = gslpp::vector<gslpp::complex>(3, 0.);
        if (flag_RI)
            RI_NLO = as_4pi_mu2 * coeffsMStoRI.transpose() * result;
        double alpha1 = as_4pi_mu2 * 4. / 3. * (12. * log(mu_2 / Mb_pole) + 6.);
        double alpha2 = as_4pi_mu2 * 4. / 3. * (6. * log(mu_2 / Mb_pole) + 13. / 2.);
        result.assign(0, -alpha2 / 2. * result(1));
        result.assign(2, -alpha1 * result(1));
        if (flag_RI)
            result += RI_NLO;
        result.assign(1, 0.);
    }
    else
    {
        throw std::runtime_error("AmpDB2::c(quark q, orders order): order not implemented");
    }
    return result;
}

◆ c_H()

gslpp::vector< gslpp::complex > AmpDB2::c_H	(	quark	q,
		orders	order
	)

private

Definition at line 1188 of file AmpDB2.cpp.

{
    gslpp::vector<gslpp::complex> result = gslpp::vector<gslpp::complex>(3, 0.);
    gslpp::complex lambda_c, lambda_u;
    switch (q)
    {
    case d:
        lambda_c = lambda_c_d;
        lambda_u = lambda_u_d;
        break;
    case s:
        lambda_c = lambda_c_s;
        lambda_u = lambda_u_s;
        break;
    default:
        throw std::runtime_error("AmpDB2::c_H(quark q): invalid quark index: ");
    }
    if (flag_RI and order == NLO)
    {
        // LO transformed to the traditional basis
        result.assign(0, -lambda_c * lambda_c * H(cc, LO) - 2. * lambda_c * lambda_u * H(cu, LO) - lambda_u * lambda_u * H(uu, LO));
        result.assign(2, -lambda_c * lambda_c * H_s(cc, LO) - 2. * lambda_c * lambda_u * H_s(cu, LO) - lambda_u * lambda_u * H_s(uu, LO));
        result.assign(0, result(0) - 0.5 * result(2));
        result.assign(1, -result(2));
        result.assign(2, 0.);
        // only NLO change to RI
        result = as_4pi_mu2 * coeffsMStoRI.transpose() * result;
        // return to tradBasis + NLO
        result.assign(0, result(0) - 0.5 * result(1) +
                             -lambda_c * lambda_c * H(cc, NLO) - 2. * lambda_c * lambda_u * H(cu, NLO) - lambda_u * lambda_u * H(uu, NLO));
        result.assign(2, -result(1) +
                             -lambda_c * lambda_c * H_s(cc, NLO) - 2. * lambda_c * lambda_u * H_s(cu, NLO) - lambda_u * lambda_u * H_s(uu, NLO));
        result.assign(1, 0.);
        return result;
    }
    if (flag_RI and order != LO)
        throw std::runtime_error("AmpDB2::c_H(quark q, orders order) order for RI not implemented");
 
    result.assign(0, -lambda_c * lambda_c * H(cc, order) - 2. * lambda_c * lambda_u * H(cu, order) - lambda_u * lambda_u * H(uu, order));
    result.assign(2, -lambda_c * lambda_c * H_s(cc, order) - 2. * lambda_c * lambda_u * H_s(cu, order) - lambda_u * lambda_u * H_s(uu, order));
    return result;
}

◆ c_H_partial()

gslpp::vector< gslpp::complex > AmpDB2::c_H_partial	(	quark	q,
		int	i
	)

private

Definition at line 1265 of file AmpDB2.cpp.

{
    gslpp::complex lambda_c, lambda_u;
    switch (q)
    {
    case d:
        lambda_c = lambda_c_d;
        lambda_u = lambda_u_d;
        break;
    case s:
        lambda_c = lambda_c_s;
        lambda_u = lambda_u_s;
        break;
    default:
        throw std::runtime_error("AmpDB2::c_H_partial(quark q, int i): invalid quark index: ");
    }
    gslpp::vector<gslpp::complex> zeros(12, 0.);
    switch (i)
    {
    case 0:
        return -lambda_c * lambda_c * H_allpartial(cc) - 2. * lambda_c * lambda_u * H_allpartial(cu) - lambda_u * lambda_u * H_allpartial(uu);
    case 1:
        return zeros;
    case 2:
        return -lambda_c * lambda_c * H_s_allpartial(cc) - 2. * lambda_c * lambda_u * H_s_allpartial(cu) - lambda_u * lambda_u * H_s_allpartial(uu);
    default:
        throw std::runtime_error("AmpDB2::c_H_partial(quark q, int i): invalid index i");
    }
}

◆ compute_deltas_1overm()

void AmpDB2::compute_deltas_1overm ( quark q )

private

Definition at line 1006 of file AmpDB2.cpp.

{
    // hep-ph/0612167 eq.24
    compute_g();
    for (quarks qq = cc; qq <= uu; qq = quarks(qq + 1))
    {
        cache_deltas_1overm[index_deltas(qq, q)] = 0.;
        for (int i = 0; i <= 3; i++)
        {
            cache_deltas_1overm[index_deltas(qq, q)] += g(qq, i) * me_R(i);
        }
        for (int i = 0; i <= 2; i++)
        {
            cache_deltas_1overm[index_deltas(qq, q)] += gtilde(qq, i) * me_Rtilde(i);
        }
    }
    return;
}

◆ compute_deltas_1overm_tradBasis()

void AmpDB2::compute_deltas_1overm_tradBasis ( quark q )

private

Definition at line 960 of file AmpDB2.cpp.

{
    // hep-ph/0308029v2 eq.20
    cache_deltas_1overm_tradBasis[index_deltas(cc, q)] = sqrt1minus4z_1overm * ((1 + 2. * z_1overm) * (K_2 * (me_R(2) + 2. * me_R(4)) - 2. * K_1 * (me_R(1) + me_R(2))) - 12. * z_1overm2 / (1. - 4. * z_1overm) * (K_1 * (me_R(2) + 2. * me_R(3)) + 2. * K_2 * me_R(3)));
    cache_deltas_1overm_tradBasis[index_deltas(cu, q)] = oneminusz_1overm2 * ((1. + 2. * z_1overm) * (K_2 * (me_R(2) + 2. * me_R(4)) - 2. * K_1 * (me_R(1) + me_R(2))) - 6. * z_1overm2 / (1. - z_1overm) * (K_1 * (me_R(2) + 2. * me_R(3)) + 2. * K_2 * me_R(3)));
    cache_deltas_1overm_tradBasis[index_deltas(uu, q)] = K_2 * (me_R(2) + 2. * me_R(4)) - 2. * K_1 * (me_R(1) + me_R(2));
    return;
}

◆ compute_g()

void AmpDB2::compute_g ( )

private

Definition at line 1025 of file AmpDB2.cpp.

{
    // hep-ph/0612167
    // equation (25)
    double cache_z = z_1overm;
    double cache_z2 = z_1overm2;
    double cache_oneminusz_1overm2 = oneminusz_1overm2;
    double cache_sqrt1minus4z = sqrt1minus4z_1overm;
    for (quarks qq = cc; qq <= uu; qq = quarks(qq + 2))
    {
        cacheg[indexg(qq, 0)] = sqrt1minus4z_1overm * (1. + 2. * z_1overm) * K_1;
        cacheg[indexg(qq, 1)] = -2. * sqrt1minus4z_1overm * (1. + 2. * z_1overm) * K_1;
        cacheg[indexg(qq, 2)] = -2. * (1. - 2. * z_1overm - 2. * z_1overm2) / sqrt1minus4z_1overm * K_1;
        cacheg[indexg(qq, 3)] = -24. * z_1overm2 / sqrt1minus4z_1overm * K_1;
        cachegtilde[indexg(qq, 0)] = -2. * sqrt1minus4z_1overm * (1. + 2. * z_1overm) * K_2;
        cachegtilde[indexg(qq, 1)] = -2. * (1. - 2. * z_1overm - 2. * z_1overm2) / sqrt1minus4z_1overm * K_2;
        cachegtilde[indexg(qq, 2)] = -24. * z_1overm2 / sqrt1minus4z_1overm * K_2;
        z_1overm = 0.;
        z_1overm2 = 0.;
        oneminusz_1overm2 = 1.;
        sqrt1minus4z_1overm = 1.;
    }
    // equation (26)
    z_1overm = cache_z;
    z_1overm2 = cache_z2;
    oneminusz_1overm2 = cache_oneminusz_1overm2;
    sqrt1minus4z_1overm = cache_sqrt1minus4z;
    cacheg[indexg(cu, 0)] = oneminusz_1overm2 * (1. + 2. * z_1overm) * K_1;
    cacheg[indexg(cu, 1)] = -2. * oneminusz_1overm2 * (1. + 2. * z_1overm) * K_1;
    cacheg[indexg(cu, 2)] = -2. * (1. - z_1overm) * (1. + z_1overm + z_1overm2) * K_1;
    ;
    cacheg[indexg(cu, 3)] = -12. * (1. - z_1overm) * z_1overm2 * K_1;
    cachegtilde[indexg(cu, 0)] = -2. * oneminusz_1overm2 * (1. + 2. * z_1overm) * K_2;
    cachegtilde[indexg(cu, 1)] = -2. * (1. - z_1overm) * (1. + z_1overm + z_1overm2) * K_2;
    cachegtilde[indexg(cu, 2)] = -12. * (1. - z_1overm) * z_1overm2 * K_2;
    return;
}

◆ compute_matrixelements()

void AmpDB2::compute_matrixelements	(	quark	q,
		orders	order
	)

private

A method to compute all DB=2 Wilson coefficients (me, me_R, me_Rtilde)

Parameters

[in]	quark	index of the neutral B mesons
[in]	order	QCD orders used for R_0 @detail requires computeCKMandMasses() before use

Definition at line 774 of file AmpDB2.cpp.

{
    double Mq;
    double Mb_mu2 = mySM.Mrun(mu_2,
                              mySM.getQuarks(QCD::BOTTOM).getMass_scale(),
                              mySM.getQuarks(QCD::BOTTOM).getMass(), QCD::BOTTOM, FULLNNLO);
    double MBq2;
    double KBq;
    double FBq2;
    switch (q)
    {
    case d:
        me = mySM.getBBd().getBpars();
        Mq = mySM.Mrun(mu_2,
                       mySM.getQuarks(QCD::DOWN).getMass_scale(),
                       mySM.getQuarks(QCD::DOWN).getMass(), QCD::DOWN, FULLNNLO);
        MBq2 = MB * MB;
        FBq2 = mySM.getMesons(QCD::B_D).getDecayconst() * mySM.getMesons(QCD::B_D).getDecayconst();
        break;
    case s:
        me = mySM.getBBs().getBpars();
        Mq = mySM.Mrun(mu_2,
                       mySM.getQuarks(QCD::STRANGE).getMass_scale(),
                       mySM.getQuarks(QCD::STRANGE).getMass(), QCD::STRANGE, FULLNNLO);
        MBq2 = MB_s * MB_s;
        FBq2 = mySM.getMesons(QCD::B_S).getDecayconst() * mySM.getMesons(QCD::B_S).getDecayconst();
        break;
    default:
        throw std::runtime_error("AmpDB2::compute_matrixelements(quark q): invalid quark index: ");
    }
    // pole mass for mbpow like in: hep-ph/0612167 eq. 28
    // double Mbpow = mySM.Mbar2Mp(mySM.getQuarks(QCD::BOTTOM).getMass());
    // double Mbpow2 = Mbpow * Mbpow;
    // KBq = MBq2 / ((Mbpow + Mq) * (Mbpow + Mq));
 
    // arXiv:1907.01025v2 equation (4)
    KBq = MBq2 / ((Mb_mu2 + Mq) * (Mb_mu2 + Mq));
    me(0) *= 8. / 3. * MBq2 * FBq2;
    me(1) *= -5. / 3. * KBq * MBq2 * FBq2;
    me(2) *= 1. / 3. * KBq * MBq2 * FBq2;
    me(3) *= 2. * (KBq + 1. / 6.) * MBq2 * FBq2;
    me(4) *= 2. / 3. * (KBq + 3. / 2.) * MBq2 * FBq2;
 
    // old parameterization from hep-ph/0612167 eq. 28 (or hep-ph/0308029 eq.26)
    /*
    me_R(0) = me(1) + me(2) + 0.5 * me(0);
    me_R(1) = Mq/Mb * me(3);
    me_R(2) = -2./3. * FBq2 * MBq2 * (MBq2 / Mbpow2 - 1.);
    me_R(3) =  7./6. * FBq2 * MBq2 * (MBq2 / Mbpow2 - 1.);
    me_R(4) = 0.5 * (me(2) + 0.5 * me(0) + me(1) - 2. * Mq/Mbpow * me(4) + me_R(2));
    */
 
    // Gerlach thesis eq.7.5, 7.6
    switch (q)
    {
    case d:
        me_R(2) = mySM.getBBd_subleading().getBpars()(0);
        me_R(3) = mySM.getBBd_subleading().getBpars()(1);
        break;
    case s:
        me_R(2) = mySM.getBBs_subleading().getBpars()(0);
        me_R(3) = mySM.getBBs_subleading().getBpars()(1);
        break;
    }
    me_R(2) *= MBq2 * FBq2;
    me_R(3) *= MBq2 * FBq2;
    me_R(4) = 0.5 * (me(2) + 0.5 * me(0) + me(1) - 2. * Mq / Mb_mu2 * me(4) + me_R(2));
 
    me_Rtilde(1) = -me_R(2);
    me_Rtilde(2) = me_R(3) + 0.5 * me_R(2);
 
    // subleading matrix elements that are not independent
    // Gerlach thesis eq. (3.64)
    double L_2 = 2. * log(mu_2 / Mb_mu2);
    double as1_me0 = 0., as1_me2 = 0.;
    if (order == FULLNLO)
    {
        as1_me0 = 4. * L_2 + 26. / 3.;
        as1_me2 = 8. * L_2 + 8.;
    }
    me_R(0) = 0.5 * (1. + as1_me0 * as_4pi_mu2) * me(0) + me(1) + (1. + as1_me2 * as_4pi_mu2) * me(2);
    // arxiv/1907.01025 eq.26
    me_R(1) = Mq / Mb_mu2 * me(3);
    me_Rtilde(0) = Mq / Mb_mu2 * me(4);
 
    // switch matrix elements to RI scheme
    if (flag_RI)
    {
        me_R(0) = 0.5 * me(0) + me(1) + me(2);
        if (order == FULLNLO)
            me += -as_4pi_mu2 * meMStoRI * me;
    }
 
    // matrix elements needed for the leading term of Gamma21overM21
    meoverme0(0) = 1.;
    meoverme0(1) = me(1) / me(0);
    meoverme0(2) = me(2) / me(0);
    return;
}

◆ compute_partialNNLO()

void AmpDB2::compute_partialNNLO ( )

private

Definition at line 2562 of file AmpDB2.cpp.

{
    for (quarks qq = cc; qq <= uu; qq = quarks(qq + 1))
    {
        for (int i = 1; i <= 6; i++)
        {
            // not all terms used for n=2
            for (int j = i; j <= 8; j++)
            {
                if (j == 3)
                    j = 8;
                if (i >= 3)
                    j = 8;
                cache_p[index_p(qq, i, j, 2)] = 0.;
                cache_ps[index_p(qq, i, j, 2)] = 0.;
            }
        }
    }
    for (int i = 0; i < 8; i++)
    {
        if (i == 6)
            i = 7;
        C_Misiak_NNLO[i] = 0.;
    }
    return;
}

◆ compute_pp_s()

void AmpDB2::compute_pp_s ( )

private

Definition at line 1457 of file AmpDB2.cpp.

{
    // input didn't change nothing to compute
    // double currentInput_compute_pp_s[4] = {z, mu_1, mu_2, mu_b}; // Removed unused variable
    // if (lastInput_compute_pp_s == currentInput_compute_pp_s) return;
 
    double cache_logz = logz;
    if (!flag_resumz)
    {
        throw std::runtime_error("AmpDB2::compute_pp_s(): only implemented for resummed z");
    }
 
    // remember value of z after setting it to 0 for calculation of uu coefficients
    const double cache_z = z;
    const double cache_sqrt1minus4z = sqrt1minus4z;
    for (quarks qq = cc; qq <= uu; qq = quarks(qq + 2))
    {
        // Gerlach thesis eq. (6.6)
        cache_p[index_p(qq, 1, 1, 0)] = (23. * sqrt1minus4z) / 72. - (43. * sqrt1minus4z * z) / 36.;
        cache_p[index_p(qq, 1, 2, 0)] = sqrt1minus4z / 6. - (5. * sqrt1minus4z * z) / 3.;
        cache_p[index_p(qq, 2, 2, 0)] = sqrt1minus4z - sqrt1minus4z * z;
        cache_ps[index_p(qq, 1, 1, 0)] = (-5. * sqrt1minus4z) / 9. - (10. * sqrt1minus4z * z) / 9.;
        cache_ps[index_p(qq, 1, 2, 0)] = (-4. * sqrt1minus4z) / 3. - (8. * sqrt1minus4z * z) / 3.;
        cache_ps[index_p(qq, 2, 2, 0)] = sqrt1minus4z + 2. * sqrt1minus4z * z;
 
        cache_p[index_p(qq, 1, 3, 0)] = 4. / 3.;
        cache_p[index_p(qq, 1, 4, 0)] = -5. / 36.;
        cache_p[index_p(qq, 1, 5, 0)] = (64. * sqrt1minus4z) / 3. - (160. * sqrt1minus4z * z) / 3.;
        cache_p[index_p(qq, 1, 6, 0)] = (-20. * sqrt1minus4z) / 9. - (4. * sqrt1minus4z * z) / 9.;
        cache_p[index_p(qq, 2, 3, 0)] = 1.;
        cache_p[index_p(qq, 2, 4, 0)] = 5. / 6.;
        cache_p[index_p(qq, 2, 5, 0)] = 16. * sqrt1minus4z - 40. * sqrt1minus4z * z;
        cache_p[index_p(qq, 2, 6, 0)] = (40. * sqrt1minus4z) / 3. + (8. * sqrt1minus4z * z) / 3.;
 
        cache_ps[index_p(qq, 1, 3, 0)] = -8. / 3.;
        cache_ps[index_p(qq, 1, 4, 0)] = -2. / 9.;
        cache_ps[index_p(qq, 1, 5, 0)] = -128. / 3.;
        cache_ps[index_p(qq, 1, 6, 0)] = -32. / 9.;
        cache_ps[index_p(qq, 2, 3, 0)] = -2.;
        cache_ps[index_p(qq, 2, 4, 0)] = 4. / 3.;
        cache_ps[index_p(qq, 2, 5, 0)] = -32.;
        cache_ps[index_p(qq, 2, 6, 0)] = 64. / 3.;
 
        // Gerlach thesis eq. (6.9)
        z = 0.;
        sqrt1minus4z = 1.;
    }
    z = cache_z;
    sqrt1minus4z = cache_sqrt1minus4z;
 
    // Gerlach thesis eq. (6.8)
    double n_l = 3.;
    double n_v = 1.;
    cache_p[index_p(cc, 3, 3, 0)] = 3. * (n_l + n_v) + 2.;
    cache_p[index_p(cc, 3, 4, 0)] = 7. / 3.;
    cache_p[index_p(cc, 3, 5, 0)] = 60. * (n_l + n_v) + 64.;
    cache_p[index_p(cc, 3, 6, 0)] = 112. / 3.;
    cache_p[index_p(cc, 4, 4, 0)] = 5. * (n_l + n_v) / 12. + 13. / 72.;
    cache_p[index_p(cc, 4, 5, 0)] = 112. / 3.;
    cache_p[index_p(cc, 4, 6, 0)] = 25. / 3. * (n_l + n_v) + 52. / 9.;
    cache_p[index_p(cc, 5, 5, 0)] = 512. + 408. * n_l + 408. * n_v * sqrt1minus4z - 480. * n_v * sqrt1minus4z * z;
    cache_p[index_p(cc, 5, 6, 0)] = 1792. / 3.;
    cache_p[index_p(cc, 6, 6, 0)] = 416. / 9. + (170. * n_l) / 3. + (170. * n_v * sqrt1minus4z) / 3. + (124. * n_v * sqrt1minus4z * z) / 3.;
 
    cache_ps[index_p(cc, 3, 3, 0)] = -6. * (n_l + n_v) - 1.;
    cache_ps[index_p(cc, 3, 4, 0)] = -8. / 3.;
    cache_ps[index_p(cc, 3, 5, 0)] = -120. * (n_l + n_v) - 32.;
    cache_ps[index_p(cc, 3, 6, 0)] = -128. / 3.;
    cache_ps[index_p(cc, 4, 4, 0)] = 2. / 3. * (n_l + n_v) - 7. / 9.;
    cache_ps[index_p(cc, 4, 5, 0)] = -128. / 3.;
    cache_ps[index_p(cc, 4, 6, 0)] = 40. / 3. * (n_l + n_v) - 224. / 9.;
    cache_ps[index_p(cc, 5, 5, 0)] = -816. * (n_l + n_v) - 256.;
    cache_ps[index_p(cc, 5, 6, 0)] = -2048. / 3.;
    cache_ps[index_p(cc, 6, 6, 0)] = 272. / 3. * (n_l + n_v) - 1792. / 9.;
 
    // Gerlach thesis eq. (6.9)
    for (int i = 3; i <= 6; i++)
    {
        for (int j = i; j <= 6; j++)
        {
            cache_p[index_p(uu, i, j, 0)] = cache_p[index_p(cc, i, j, 0)];
            cache_ps[index_p(uu, i, j, 0)] = cache_ps[index_p(cc, i, j, 0)];
        }
    }
 
    // Gerlach thesis eq. (6.10)
    for (int i = 1; i <= 6; i++)
    {
        for (int j = i; j <= 6; j++)
        {
            cache_p[index_p(cu, i, j, 0)] =
                0.5 * (cache_p[index_p(cc, i, j, 0)] + cache_p[index_p(uu, i, j, 0)]);
            cache_ps[index_p(cu, i, j, 0)] =
                0.5 * (cache_ps[index_p(cc, i, j, 0)] + cache_ps[index_p(uu, i, j, 0)]);
        }
    }
    double L_1 = 2. * log(mu_1 / Mb);
    double L_2 = 2. * log(mu_2 / Mb);
    double log_some_z = log((2 * sqrt1minus4z) / (1 + sqrt1minus4z));
 
    cache_p[index_p(cc, 1, 1, 1)] = (169. * Dilogsigma) / 27. + (92. * Dilogsigma2) / 27. + (92. * log2sigma) / 27. - (191. * logsigma) / 54. +
                                    (46. * log1minus4z * logsigma) / 27. + (169. * logsigma * log_some_z) / 54. - (92. * logsigma * logz) / 27. +
                                    (1211. * sqrt1minus4z) / 324. + (19. * L_1 * sqrt1minus4z) / 18. + (149. * L_2 * sqrt1minus4z) / 108. -
                                    (8. * log1minus4z * sqrt1minus4z) / 3. + (43. * logz * sqrt1minus4z) / 12. - (5. * logsigma) / (216. * z) +
                                    (5. * logz * sqrt1minus4z) / (216. * z) - (340. * Dilogsigma * z) / 9. - 19. * Dilogsigma2 * z -
                                    19. * log2sigma * z + (371. * logsigma * z) / 27. - (19. * log1minus4z * logsigma * z) / 2. -
                                    (170. * logsigma * log_some_z * z) / 9. + 19. * logsigma * logz * z - (6917. * sqrt1minus4z * z) / 324. -
                                    (23. * L_1 * sqrt1minus4z * z) / 9. - (49. * L_2 * sqrt1minus4z * z) / 54. +
                                    (128. * log1minus4z * sqrt1minus4z * z) / 9. - (1951. * logz * sqrt1minus4z * z) / 108. +
                                    (1364. * Dilogsigma * z2) / 27. + (682. * Dilogsigma2 * z2) / 27. + (682. * log2sigma * z2) / 27. -
                                    (263. * logsigma * z2) / 54. + (341. * log1minus4z * logsigma * z2) / 27. +
                                    (682. * logsigma * log_some_z * z2) / 27. - (682. * logsigma * logz * z2) / 27. -
                                    (295. * sqrt1minus4z * z2) / 54. + (295. * logsigma * z3) / 27. - (5. * (M_PI2)) / 108. + (z * (M_PI2)) / 54.;
    cache_p[index_p(cc, 1, 2, 1)] = (92. * Dilogsigma) / 9. + (16. * Dilogsigma2) / 9. + (16. * log2sigma) / 9. - (275. * logsigma) / 36. +
                                    (8. * log1minus4z * logsigma) / 9. + (46. * logsigma * log_some_z) / 9. - (16. * logsigma * logz) / 9. -
                                    (476. * sqrt1minus4z) / 27. - (37. * L_1 * sqrt1minus4z) / 6. + (19. * L_2 * sqrt1minus4z) / 9. +
                                    (27. * logz * sqrt1minus4z) / 4. + (4. * logsigma) / (9. * z) - (4. * logz * sqrt1minus4z) / (9. * z) -
                                    (176. * Dilogsigma * z) / 3. - 28. * Dilogsigma2 * z - 28. * log2sigma * z + (815. * logsigma * z) / 18. -
                                    14. * log1minus4z * logsigma * z - (88. * logsigma * log_some_z * z) / 3. + 28. * logsigma * logz * z +
                                    (917. * sqrt1minus4z * z) / 27. + (41. * L_1 * sqrt1minus4z * z) / 3. + (2. * L_2 * sqrt1minus4z * z) / 9. +
                                    (64. * log1minus4z * sqrt1minus4z * z) / 3. - (392. * logz * sqrt1minus4z * z) / 9. + (688. * Dilogsigma * z2) / 9. +
                                    (344. * Dilogsigma2 * z2) / 9. + (344. * log2sigma * z2) / 9. - (269. * logsigma * z2) / 9. +
                                    (172. * log1minus4z * logsigma * z2) / 9. + (344. * logsigma * log_some_z * z2) / 9. - (344. * logsigma * logz * z2) / 9. -
                                    (91. * sqrt1minus4z * z2) / 9. + (182. * logsigma * z3) / 9. + (5. * (M_PI2)) / 9. - (2. * z * (M_PI2)) / 9.;
    cache_p[index_p(cc, 2, 2, 1)] = (4. * Dilogsigma) / 3. + (32. * Dilogsigma2) / 3. + (32. * log2sigma) / 3. - (41. * logsigma) / 6. +
                                    (16. * log1minus4z * logsigma) / 3. + (2. * logsigma * log_some_z) / 3. - (32. * logsigma * logz) / 3. +
                                    (103. * sqrt1minus4z) / 18. - L_1 * sqrt1minus4z + (2. * L_2 * sqrt1minus4z) / 3. -
                                    12. * log1minus4z * sqrt1minus4z + (21. * logz * sqrt1minus4z) / 2. - (11. * logsigma) / (6. * z) +
                                    (11. * logz * sqrt1minus4z) / (6. * z) - 16. * Dilogsigma * z - 12. * Dilogsigma2 * z - 12. * log2sigma * z +
                                    (77. * logsigma * z) / 3. - 6. * log1minus4z * logsigma * z - 8. * logsigma * log_some_z * z + 12. * logsigma * logz * z +
                                    (34. * sqrt1minus4z * z) / 9. + 10. * L_1 * sqrt1minus4z * z - (14. * L_2 * sqrt1minus4z * z) / 3. +
                                    8. * log1minus4z * sqrt1minus4z * z - (73. * logz * sqrt1minus4z * z) / 3. + (80. * Dilogsigma * z2) / 3. +
                                    (40. * Dilogsigma2 * z2) / 3. + (40. * log2sigma * z2) / 3. - (16. * logsigma * z2) / 3. +
                                    (20. * log1minus4z * logsigma * z2) / 3. + (40. * logsigma * log_some_z * z2) / 3. - (40. * logsigma * logz * z2) / 3. +
                                    (16. * sqrt1minus4z * z2) / 3. - (32. * logsigma * z3) / 3. - (5. * (M_PI2)) / 3. + (2. * z * (M_PI2)) / 3.;
    cache_ps[index_p(cc, 1, 1, 1)] = (-344. * Dilogsigma) / 27. - (160. * Dilogsigma2) / 27. - (160. * log2sigma) / 27. + (320. * logsigma) / 27. -
                                     (80. * log1minus4z * logsigma) / 27. - (172. * logsigma * log_some_z) / 27. + (160. * logsigma * logz) / 27. -
                                     (10. * sqrt1minus4z) / 81. - (4. * L_1 * sqrt1minus4z) / 9. - (40. * L_2 * sqrt1minus4z) / 27. +
                                     (80. * log1minus4z * sqrt1minus4z) / 9. - 12. * logz * sqrt1minus4z + (2. * logsigma) / (27. * z) -
                                     (2. * logz * sqrt1minus4z) / (27. * z) + (8. * Dilogsigma2 * z) / 9. + (8. * log2sigma * z) / 9. -
                                     (214. * logsigma * z) / 27. + (4. * log1minus4z * logsigma * z) / 9. - (8. * logsigma * logz * z) / 9. -
                                     (1511. * sqrt1minus4z * z) / 81. - (8. * L_1 * sqrt1minus4z * z) / 9. - (80. * L_2 * sqrt1minus4z * z) / 27. +
                                     (160. * log1minus4z * sqrt1minus4z * z) / 9. - (610. * logz * sqrt1minus4z * z) / 27. +
                                     (1376. * Dilogsigma * z2) / 27. + (688. * Dilogsigma2 * z2) / 27. + (688. * log2sigma * z2) / 27. -
                                     (460. * logsigma * z2) / 27. + (344. * log1minus4z * logsigma * z2) / 27. + (688. * logsigma * log_some_z * z2) / 27. -
                                     (688. * logsigma * logz * z2) / 27. - (590. * sqrt1minus4z * z2) / 27. + (1180. * logsigma * z3) / 27. -
                                     (2. * (M_PI2)) / 27. - (4. * z * (M_PI2)) / 27.;
    cache_ps[index_p(cc, 1, 2, 1)] = (-160. * Dilogsigma) / 9. - (128. * Dilogsigma2) / 9. - (128. * log2sigma) / 9. + (238. * logsigma) / 9. -
                                     (64. * log1minus4z * logsigma) / 9. - (80. * logsigma * log_some_z) / 9. + (128. * logsigma * logz) / 9. +
                                     (100. * sqrt1minus4z) / 27. + (4. * L_1 * sqrt1minus4z) / 3. - (32. * L_2 * sqrt1minus4z) / 9. +
                                     (64. * log1minus4z * sqrt1minus4z) / 3. - 26. * logz * sqrt1minus4z - (2. * logsigma) / (9. * z) +
                                     (2. * logz * sqrt1minus4z) / (9. * z) - (32. * Dilogsigma2 * z) / 3. - (32. * log2sigma * z) / 3. +
                                     (16. * logsigma * z) / 9. - (16. * log1minus4z * logsigma * z) / 3. + (32. * logsigma * logz * z) / 3. -
                                     (442. * sqrt1minus4z * z) / 27. + (8. * L_1 * sqrt1minus4z * z) / 3. - (64. * L_2 * sqrt1minus4z * z) / 9. +
                                     (128. * log1minus4z * sqrt1minus4z * z) / 3. - (560. * logz * sqrt1minus4z * z) / 9. + (640. * Dilogsigma * z2) / 9. +
                                     (320. * Dilogsigma2 * z2) / 9. + (320. * log2sigma * z2) / 9. - (728. * logsigma * z2) / 9. +
                                     (160. * log1minus4z * logsigma * z2) / 9. + (320. * logsigma * log_some_z * z2) / 9. - (320. * logsigma * logz * z2) / 9. -
                                     (364. * sqrt1minus4z * z2) / 9. + (728. * logsigma * z3) / 9. + (8. * (M_PI2)) / 9. + (16. * z * (M_PI2)) / 9.;
    cache_ps[index_p(cc, 2, 2, 1)] = (-32. * Dilogsigma) / 3. + (32. * Dilogsigma2) / 3. + (32. * log2sigma) / 3. - (28. * logsigma) / 3. +
                                     (16. * log1minus4z * logsigma) / 3. - (16. * logsigma * log_some_z) / 3. - (32. * logsigma * logz) / 3. +
                                     (272. * sqrt1minus4z) / 9. + 8. * L_1 * sqrt1minus4z + (8. * L_2 * sqrt1minus4z) / 3. -
                                     16. * log1minus4z * sqrt1minus4z + 12. * logz * sqrt1minus4z - (4. * logsigma) / (3. * z) +
                                     (4. * logz * sqrt1minus4z) / (3. * z) + 32. * Dilogsigma2 * z + 32. * log2sigma * z - (40. * logsigma * z) / 3. +
                                     16. * log1minus4z * logsigma * z - 32. * logsigma * logz * z + (664. * sqrt1minus4z * z) / 9. +
                                     16. * L_1 * sqrt1minus4z * z + (16. * L_2 * sqrt1minus4z * z) / 3. - 32. * log1minus4z * sqrt1minus4z * z +
                                     (104. * logz * sqrt1minus4z * z) / 3. + (128. * Dilogsigma * z2) / 3. + (64. * Dilogsigma2 * z2) / 3. +
                                     (64. * log2sigma * z2) / 3. + (272. * logsigma * z2) / 3. + (32. * log1minus4z * logsigma * z2) / 3. +
                                     (64. * logsigma * log_some_z * z2) / 3. - (64. * logsigma * logz * z2) / 3. + (64. * sqrt1minus4z * z2) / 3. -
                                     (128. * logsigma * z3) / 3. - (8. * (M_PI2)) / 3. - (16. * z * (M_PI2)) / 3.;
 
    cache_p[index_p(uu, 1, 1, 1)] = 299. / 81. + (19. * L_1) / 18. + (149. * L_2) / 108. - (5. * (M_PI2)) / 108.;
    cache_p[index_p(uu, 1, 2, 1)] = -452. / 27. - (37. * L_1) / 6. + (19. * L_2) / 9. + (5. * (M_PI2)) / 9.;
    cache_p[index_p(uu, 2, 2, 1)] = 37. / 18. - L_1 + (2. * L_2) / 3. - (5. * (M_PI2)) / 3.;
    cache_ps[index_p(uu, 1, 1, 1)] = 2. / 81. - (4. * L_1) / 9. - (40. * L_2) / 27. - (2. * (M_PI2)) / 27.;
    cache_ps[index_p(uu, 1, 2, 1)] = 88. / 27. + (4. * L_1) / 3. - (32. * L_2) / 9. + (8. * (M_PI2)) / 9.;
    cache_ps[index_p(uu, 2, 2, 1)] = 248. / 9. + 8. * L_1 + (8. * L_2) / 3. - (8. * (M_PI2)) / 3.;
 
    // pengFlag terms from P12-36 here in full z dependence
    cache_p[index_p(cc, 1, 1, 1)] += (-5. * logsigma) / 324. - (5. * sqrt1minus4z) / 486. -
                                     (5. * L_1 * sqrt1minus4z) / 324. + (5. * logz * sqrt1minus4z) / 324. -
                                     (20. * sqrt1minus4z * z) / 243. - (5. * L_1 * sqrt1minus4z * z) / 162. +
                                     (5. * logz * sqrt1minus4z * z) / 162. + (5. * logsigma * z2) / 27. -
                                     (10. * sqrt1minus4z * z2) / 81. + (20. * logsigma * z3) / 81.;
    cache_p[index_p(cc, 1, 2, 1)] += (5. * logsigma) / 27. + (10. * sqrt1minus4z) / 81. + (5. * L_1 * sqrt1minus4z) / 27. -
                                     (5. * logz * sqrt1minus4z) / 27. + (80. * sqrt1minus4z * z) / 81. +
                                     (10. * L_1 * sqrt1minus4z * z) / 27. - (10. * logz * sqrt1minus4z * z) / 27. -
                                     (20. * logsigma * z2) / 9. + (40. * sqrt1minus4z * z2) / 27. -
                                     (80. * logsigma * z3) / 27.;
    cache_p[index_p(cc, 2, 2, 1)] += (-5. * logsigma) / 9. - (10. * sqrt1minus4z) / 27. - (5. * L_1 * sqrt1minus4z) / 9. +
                                     (5. * logz * sqrt1minus4z) / 9. - (80. * sqrt1minus4z * z) / 27. -
                                     (10. * L_1 * sqrt1minus4z * z) / 9. + (10. * logz * sqrt1minus4z * z) / 9. +
                                     (20. * logsigma * z2) / 3. - (40. * sqrt1minus4z * z2) / 9. + (80. * logsigma * z3) / 9.;
    cache_ps[index_p(cc, 1, 1, 1)] += (-2. * logsigma) / 81. - (4. * sqrt1minus4z) / 243. -
                                      (2. * L_1 * sqrt1minus4z) / 81. + (2. * logz * sqrt1minus4z) / 81. -
                                      (32. * sqrt1minus4z * z) / 243. - (4. * L_1 * sqrt1minus4z * z) / 81. +
                                      (4. * logz * sqrt1minus4z * z) / 81. + (8. * logsigma * z2) / 27. -
                                      (16. * sqrt1minus4z * z2) / 81. + (32. * logsigma * z3) / 81.;
    cache_ps[index_p(cc, 1, 2, 1)] += (8. * logsigma) / 27. + (16. * sqrt1minus4z) / 81. + (8. * L_1 * sqrt1minus4z) / 27. -
                                      (8. * logz * sqrt1minus4z) / 27. + (128. * sqrt1minus4z * z) / 81. +
                                      (16. * L_1 * sqrt1minus4z * z) / 27. - (16. * logz * sqrt1minus4z * z) / 27. -
                                      (32. * logsigma * z2) / 9. + (64. * sqrt1minus4z * z2) / 27. -
                                      (128. * logsigma * z3) / 27.;
    cache_ps[index_p(cc, 2, 2, 1)] += (-8. * logsigma) / 9. - (16. * sqrt1minus4z) / 27. - (8. * L_1 * sqrt1minus4z) / 9. +
                                      (8. * logz * sqrt1minus4z) / 9. - (128. * sqrt1minus4z * z) / 27. -
                                      (16. * L_1 * sqrt1minus4z * z) / 9. + (16. * logz * sqrt1minus4z * z) / 9. +
                                      (32. * logsigma * z2) / 3. - (64. * sqrt1minus4z * z2) / 9. + (128. * logsigma * z3) / 9.;
 
    cache_p[index_p(uu, 1, 1, 1)] += -5. / 486. - (5. * L_1) / 324.;
    cache_p[index_p(uu, 1, 2, 1)] += 10. / 81. + (5. * L_1) / 27.;
    cache_p[index_p(uu, 2, 2, 1)] += -10. / 27. - (5. * L_1) / 9.;
    cache_ps[index_p(uu, 1, 1, 1)] += -4. / 243. - (2. * L_1) / 81.;
    cache_ps[index_p(uu, 1, 2, 1)] += 16. / 81. + (8. * L_1) / 27.;
    cache_ps[index_p(uu, 2, 2, 1)] += -16. / 27. - (8. * L_1) / 9.;
 
    // Corrections from LO coefficients when resummation logz
    double z_replace = -6. * 4. / 3. * logz * z;
    cache_p[index_p(cc, 1, 1, 1)] += (-11. / 6. + 43. / 6. * z) / sqrt1minus4z * z_replace;
    cache_p[index_p(cc, 1, 2, 1)] += (-2. + 10. * z) / sqrt1minus4z * z_replace;
    cache_p[index_p(cc, 2, 2, 1)] += (-3. + 6. * z) / sqrt1minus4z * z_replace;
    cache_ps[index_p(cc, 1, 1, 1)] += (20. / 3. * z) / sqrt1minus4z * z_replace;
    cache_ps[index_p(cc, 1, 2, 1)] += (16. * z) / sqrt1minus4z * z_replace;
    cache_ps[index_p(cc, 2, 2, 1)] += (-12. * z) / sqrt1minus4z * z_replace;
 
    for (quarks qq = cc; qq <= uu; qq = quarks(qq + 2))
    {
        // equation (6.14)
        cache_p[index_p(qq, 1, 3, 1)] = (320. / 9. - 4. * L_1) * z + 47. / 18. * L_1 + 56. / 9. * L_2 - 5. / (18. * sqrt3) * M_PI + 1523. / 108.;
        cache_p[index_p(qq, 1, 4, 1)] = (59. / 3. * L_1 + 5. / 9. * M_PI2 + 4565. / 108.) * z - 281. / 108. * L_1 + L_2 / 54. + 5. / 18. * M_PI2 - 25. / (108. * sqrt3) * M_PI - 712. / 81.;
        cache_p[index_p(qq, 1, 5, 1)] = z * (-136. * L_1 - 192. * L_2 - 768. * logz - 16408. / 9.) + 376. / 9. * L_1 + 896. / 9. * L_2 - 40. / (9. * sqrt3) * M_PI + 318.;
        cache_p[index_p(qq, 1, 6, 1)] = z * (764. / 3. * L_1 + 8. * L_2 + 32. * logz + 8. / 9. * M_PI2 + 22850. / 27.) - 1259. / 27. * L_1 + 8. / 27. * L_2 + 40. / 9. * M_PI2 - 55. / (27. * sqrt3) * M_PI - 4243. / 27.;
        cache_p[index_p(qq, 2, 3, 1)] = (24. * L_1 + 170. / 3.) * z - 47. / 3. * L_1 + 14. / 3. * L_2 + 5. / (3. * sqrt3) * M_PI - 677. / 18.;
        cache_p[index_p(qq, 2, 4, 1)] = (26. * L_1 - 10. / 3. * M_PI2 + 1429. / 18.) * z - 35. / 9. * L_1 - L_2 / 9. - 5. / 3. * M_PI2 + 25. / (18. * sqrt3) * M_PI - 88. / 27.;
        cache_p[index_p(qq, 2, 5, 1)] = z * (816. * L_1 - 144. * L_2 - 576. * logz + 3656. / 3.) - 752. / 3. * L_1 + 224. / 3. * L_2 + 80. / (3. * sqrt3) * M_PI - 580.;
        cache_p[index_p(qq, 2, 6, 1)] = z * (128. * L_1 - 48. * L_2 - 192. * logz - 16. / 3. * M_PI2 + 6140. / 9.) - 290. / 9. * L_1 - 16. / 9. * L_2 - 80. / 3. * M_PI2 + 110. / (9. * sqrt3) * M_PI - 442. / 9.;
 
        cache_ps[index_p(qq, 1, 3, 1)] = -4. / 3. * L_1 - 64. / 9. * L_2 - 1720. / 9. * z - 4. / (9. * sqrt3) * M_PI - 130. / 27.;
        cache_ps[index_p(qq, 1, 4, 1)] = 2. * L_1 - 16. / 27. * L_2 + (80. / 27. + 8. / 9. * M_PI2) * z + 4. / 9. * M_PI2 - 10. / (27. * sqrt3) * M_PI + 404. / 81.;
        cache_ps[index_p(qq, 1, 5, 1)] = -64. / 3. * L_1 - 1024. / 9. * L_2 - 27952. / 9. * z - 64. / (9. * sqrt3) * M_PI - 2128. / 9.;
        cache_ps[index_p(qq, 1, 6, 1)] = 24. * L_1 - 256. / 27. * L_2 + (128. / 9. * M_PI2 - 520. / 27.) * z - 64. / 9. * M_PI2 - 88. / (27. * sqrt3) * M_PI + 2824. / 27.;
        cache_ps[index_p(qq, 2, 3, 1)] = 8. * L_1 - 16. / 3. * L_2 - 448. / 3. * z + 8. / (3. * sqrt3) * M_PI + 116. / 9.;
        cache_ps[index_p(qq, 2, 4, 1)] = 32. / 9. * L_2 + (488. / 9. - 16. / 3. * M_PI2) * z - 8. / 3. * M_PI2 + 20. / (9. * sqrt3) * M_PI + 272. / 27.;
        cache_ps[index_p(qq, 2, 5, 1)] = 128. * L_1 - 256. / 3. * L_2 - 6304. / 3. * z + 128. / (3. * sqrt3) * M_PI + 32. / 3.;
        cache_ps[index_p(qq, 2, 6, 1)] = 48. * L_1 + 512. / 9. * L_2 + (7520. / 9. - 256. / 3. * M_PI2) * z - 128. / 3. * M_PI2 + 176. / (9. * sqrt3) * M_PI + 1840. / 9.;
 
        // Gerlach thesis eq. (6.13) and (6.15)
        z = 0.;
    }
    z = cache_z;
    // Gerlach thesis eq. (6.15)
    cache_p[index_p(uu, 1, 4, 1)] += 5. / 9. * z;
    cache_p[index_p(uu, 1, 6, 1)] += 50. / 9. * z;
    cache_p[index_p(uu, 2, 4, 1)] += -10. / 3. * z;
    cache_p[index_p(uu, 2, 6, 1)] += -100. / 3. * z;
 
    cache_ps[index_p(uu, 1, 4, 1)] += 8. / 9. * z;
    cache_ps[index_p(uu, 1, 6, 1)] += 80. / 9. * z;
    cache_ps[index_p(uu, 2, 4, 1)] += -16. / 3. * z;
    cache_ps[index_p(uu, 2, 6, 1)] += -160. / 3. * z;
 
    // Corrections from LO coefficients when resummation logz
    cache_p[index_p(cc, 1, 5, 1)] += -96. * z_replace;
    cache_p[index_p(cc, 1, 6, 1)] += 4. * z_replace;
    cache_p[index_p(cc, 2, 5, 1)] += -72. * z_replace;
    cache_p[index_p(cc, 2, 6, 1)] += -24. * z_replace;
 
    // Gerlach thesis eq. (6.13) and (6.16)
    for (int i = 1; i <= 2; i++)
    {
        for (int j = i; j <= 6; j++)
        {
            cache_p[index_p(cu, i, j, 1)] =
                0.5 * (cache_p[index_p(cc, i, j, 1)] + cache_p[index_p(uu, i, j, 1)]);
            cache_ps[index_p(cu, i, j, 1)] =
                0.5 * (cache_ps[index_p(cc, i, j, 1)] + cache_ps[index_p(uu, i, j, 1)]);
        }
    }
 
    cache_p[index_p(cu, 1, 1, 1)] = (15. * log1minusz + 2407. * z + 894. * L_2 * z - 1065. * log1minusz * z + 1014. * log1minusz * logz * z - 1539. * log1minusz * z4 * z +
                                     1539. * logz * z4 * z - 36. * L_1 * oneminusz2 * z * (-19. + 4. * z) - 9255. * z2 - 1188. * L_2 * z2 + 5679. * log1minusz * z2 - 3738. * logz * z2 -
                                     2970. * log1minusz * logz * z2 + 10008. * z3 - 306. * L_2 * z3 - 9762. * log1minusz * z3 + 7794. * logz * z3 + 3060. * log1minusz * logz * z3 -
                                     12. * Dilogz * z * (-169. + 495. * z - 510. * z2 + 184. * z3) - 3160. * z4 + 600. * L_2 * z4 + 6672. * log1minusz * z4 - 6876. * logz * z4 -
                                     1104. * log1minusz * logz * z4 - 30. * z * (M_PI2) + 6. * z2 * (M_PI2) + 24. * z3 * (M_PI2)) /
                                    (648. * z);
    cache_p[index_p(cu, 1, 2, 1)] = -1. / 108. * (48. * log1minusz + 1856. * z - 228. * L_2 * z - 627. * log1minusz * z - 276. * log1minusz * logz * z + 378. * log1minusz * z4 * z - 378. * logz * z4 * z - 18. * L_1 * oneminusz2 * z * (-37. + 4. * z) - 3588. * z2 + 216. * L_2 * z2 + 72. * log1minusz * z2 + 588. * logz * z2 + 972. * log1minusz * logz * z2 + 1143. * z3 + 252. * L_2 * z3 + 1923. * log1minusz * z3 - 2421. * logz * z3 - 792. * log1minusz * logz * z3 + 24. * Dilogz * z * (-23. + 81. * z - 66. * z2 + 8. * z3) + 589. * z4 - 240. * L_2 * z4 - 1794. * log1minusz * z4 + 1746. * logz * z4 + 96. * log1minusz * logz * z4 - 60. * z * (M_PI2) + 12. * z2 * (M_PI2) + 48. * z3 * (M_PI2)) / z;
    cache_p[index_p(cu, 2, 2, 1)] = (70. + 12. * L_2 - 210. * log1minusz + 6. * log1minusz * logz + (33. * log1minusz) / z - 3. * z - 54. * L_2 * z + 225. * log1minusz * z - 210. * logz * z +
                                     54. * log1minusz * logz * z + 18. * L_1 * oneminusz2 * (-1. + 4. * z) - 360. * z2 + 72. * L_2 * z2 + 21. * log1minusz * z2 + 153. * logz * z2 +
                                     36. * log1minusz * logz * z2 + 293. * z3 - 30. * L_2 * z3 - 42. * log1minusz * z3 - 126. * logz * z3 - 96. * log1minusz * logz * z3 -
                                     12. * Dilogz * (-1. - 9. * z - 6. * z2 + 16. * z3) - 27. * log1minusz * z4 + 27. * logz * z4 - 30. * (M_PI2) + 6. * z * (M_PI2) + 24. * z2 * (M_PI2)) /
                                    18.;
    cache_ps[index_p(cu, 1, 1, 1)] = -1. / 162. * (12. * log1minusz + 8. * z + 240. * L_2 * z - 978. * log1minusz * z + 516. * log1minusz * logz * z + 774. * log1minusz * z4 * z - 774. * logz * z4 * z + 72. * L_1 * oneminusz2 * z * (1. + 2. * z) + 1461. * z2 - 1674. * log1minusz * z2 + 516. * logz * z2 - 36. * log1minusz * logz * z2 - 3654. * z3 - 720. * L_2 * z3 + 7008. * log1minusz * z3 - 5796. * logz * z3 - 1584. * log1minusz * logz * z3 + 24. * Dilogz * z * (43. - 3. * z - 132. * z2 + 92. * z3) + 2185. * z4 + 480. * L_2 * z4 - 5142. * log1minusz * z4 + 5346. * logz * z4 + 1104. * log1minusz * logz * z4 + 12. * z * (M_PI2) + 12. * z2 * (M_PI2)-24. * z3 * (M_PI2)) / z;
    cache_ps[index_p(cu, 1, 2, 1)] = (6. * log1minusz + 94. * z - 96. * L_2 * z + 402. * log1minusz * z - 120. * log1minusz * logz * z - 180. * log1minusz * z4 * z +
                                      180. * logz * z4 * z + 36. * L_1 * oneminusz2 * z * (1. + 2. * z) - 363. * z2 + 216. * log1minusz * z2 - 120. * logz * z2 -
                                      72. * log1minusz * logz * z2 + 792. * z3 + 288. * L_2 * z3 - 1842. * log1minusz * z3 + 1638. * logz * z3 +
                                      288. * log1minusz * logz * z3 - 48. * Dilogz * z * (5. + 3. * z - 12. * z2 + 4. * z3) - 523. * z4 - 192. * L_2 * z4 +
                                      1398. * log1minusz * z4 - 1350. * logz * z4 - 96. * log1minusz * logz * z4 + 24. * z * (M_PI2) + 24. * z2 * (M_PI2)-48. * z3 * (M_PI2)) /
                                     (27. * z);
    cache_ps[index_p(cu, 2, 2, 1)] = (4. * (12. * oneminusz2 * (1. + 2. * z) + 18. * L_1 * oneminusz2 * (1. + 2. * z) - 12. * L_2 * oneminusz2 * (1. + 2. * z) +
                                            18. * (1. + L_2) * oneminusz2 * (1. + 2. * z) - 6. * (2. * Dilogz + log1minusz * logz) * (-1. + z) * (1. + 2. * z) * (-1. + 4. * z) - (3. * log1minusz * oneminusz2 * (1. + z) * (-1. + 13. * z + 3. * z2)) / z +
                                            7. * (-1. + z) * (-5. - 8. * z + 19. * z2) + 3. * logz * z * (-2. + 3. * z - 18. * z2 + 3. * z3) +
                                            6. * (-1. + z) * (1. + 2. * z) * (M_PI2))) /
                                     9.;
 
    // pengFlag terms from P12-36 here in full z dependence
    cache_p[index_p(cu, 1, 1, 1)] += (5. * (-2. - 3. * L_1 + 3. * logz - 12. * z - 2. * sqrt1minus4z * (1. + 2. * z) -
                                            3. * L_1 * sqrt1minus4z * (1. + 2. * z) - 3. * logsigma * sqrt1minus4z * (1. + 2. * z))) /
                                     1944.;
    cache_p[index_p(cu, 1, 2, 1)] += (5. * (2. + 3. * L_1 - 3. * logz + 12. * z + 2. * sqrt1minus4z * (1. + 2. * z) +
                                            3. * L_1 * sqrt1minus4z * (1. + 2. * z) + 3. * logsigma * sqrt1minus4z * (1. + 2. * z))) /
                                     162.;
    cache_p[index_p(cu, 2, 2, 1)] += (5. * (-2. - 3. * L_1 + 3. * logz - 12. * z - 2. * sqrt1minus4z * (1. + 2. * z) -
                                            3. * L_1 * sqrt1minus4z * (1. + 2. * z) - 3. * logsigma * sqrt1minus4z * (1. + 2. * z))) /
                                     54.;
    cache_ps[index_p(cu, 1, 1, 1)] += (-2. - 3. * L_1 + 3. * logz - 12. * z - 2. * sqrt1minus4z * (1. + 2. * z) -
                                       3. * L_1 * sqrt1minus4z * (1. + 2. * z) - 3. * logsigma * sqrt1minus4z * (1. + 2. * z)) /
                                      243.;
    cache_ps[index_p(cu, 1, 2, 1)] += (4. * (2. + 3. * L_1 - 3. * logz + 12. * z + 2. * sqrt1minus4z * (1. + 2. * z) +
                                             3. * L_1 * sqrt1minus4z * (1. + 2. * z) + 3. * logsigma * sqrt1minus4z * (1. + 2. * z))) /
                                      81.;
    cache_ps[index_p(cu, 2, 2, 1)] += (4. * (-2. - 3. * L_1 + 3. * logz - 12. * z - 2. * sqrt1minus4z * (1. + 2. * z) -
                                             3. * L_1 * sqrt1minus4z * (1. + 2. * z) - 3. * logsigma * sqrt1minus4z * (1. + 2. * z))) /
                                      27.;
 
    // Corrections from LO coefficients when resummation logz
    cache_p[index_p(cu, 1, 1, 1)] += (-11. / 12. + 7. / 4. * z + 5. / 6. * z2) * z_replace;
    cache_p[index_p(cu, 1, 2, 1)] += (-1. + 3. * z - 2. * z2) * z_replace;
    cache_p[index_p(cu, 2, 2, 1)] += (-3. / 2. + 3. / 2. * z2) * z_replace;
    cache_ps[index_p(cu, 1, 1, 1)] += (10. / 3. * z - 10. / 3. * z2) * z_replace;
    cache_ps[index_p(cu, 1, 2, 1)] += (8. * z - 8. * z2) * z_replace;
    cache_ps[index_p(cu, 2, 2, 1)] += (-6. * z + 6. * z2) * z_replace;
 
    // Gerlach thesis eq. (6.18)
    cache_p[index_p(cc, 3, 3, 1)] = -154. / 9. * L_1 + 184. / 3. * L_2 + 90. * z - 5. / 3. * M_PI2 + 5. / (3. * sqrt3) * M_PI + 1390. / 27.;
    cache_p[index_p(cc, 3, 4, 1)] = -811. / 54. * L_1 + 74. / 9. * L_2 - 10. / 3. * z - 10. / 9. * M_PI2 + 70. / (9. * sqrt3) * M_PI - 27991. / 324.;
    cache_p[index_p(cc, 3, 5, 1)] = -4928. / 9. * L_1 + 3872. / 3. * L_2 + 1800. * z - 160. / 3. * M_PI2 + 160. / (3. * sqrt3) * M_PI + 16880. / 27.;
    cache_p[index_p(cc, 3, 6, 1)] = (144. * L_1 + 440. / 3.) * z - 12932. / 27. * L_1 + 1184. / 9. * L_2 - 160. / 9. * M_PI2 + 670. / (9. * sqrt3) * M_PI - 131410. / 81.;
    cache_p[index_p(cc, 4, 4, 1)] = 181. / 162. * L_1 + 127. / 108. * L_2 + (323. / 36. - 5. / 3. * M_PI2) * z - 335. / 108. * M_PI2 + 575. / (108. * sqrt3) * M_PI + 779. / 486.;
    cache_p[index_p(cc, 4, 5, 1)] = (576. * L_1 + 3836. / 3.) * z - 14912. / 27. * L_1 + 1184. / 9. * L_2 - 160. / 9. * M_PI2 + 1120. / (9. * sqrt3) * M_PI - 127990. / 81.;
    cache_p[index_p(cc, 4, 6, 1)] = (60. * L_1 - 100. / 3. * M_PI2 + 2455. / 9.) * z - 8759. / 81. * L_1 + 1088. / 27. * L_2 - 1600. / 27. * M_PI2 + 2665. / (27. * sqrt3) * M_PI - 50083. / 243.;
    cache_p[index_p(cc, 5, 5, 1)] = z * (-2592. * L_2 - 10368. * logz - 33120.) - 39424. / 9. * L_1 + 26944. / 3. * L_2 - 1280. / 3. * M_PI2 + 1280. / (3. * sqrt3) * M_PI + 347104. / 27.;
    cache_p[index_p(cc, 5, 6, 1)] = (7200. * L_1 + 74000. / 3.) * z - 240608. / 27. * L_1 + 18944. / 9. * L_2 - 2560. / 9. * M_PI2 + 10720. / (9. * sqrt3) * M_PI - 2253568. / 81.;
    cache_p[index_p(cc, 6, 6, 1)] = z * (-48. * L_1 - 144. * L_2 - 576. * logz - 248. / 3. * M_PI2 + 12290. / 9.) - 59632. / 81. * L_1 + 8848. / 27. * L_2 - 10640. / 27. * M_PI2 + 12320. / (27. * sqrt3) * M_PI - 662144. / 243.;
 
    cache_ps[index_p(cc, 3, 3, 1)] = 176. / 9. * L_1 - 200. / 3. * L_2 - 432. * z - 8. / 3. * M_PI2 + 8. / (3. * sqrt3) * M_PI - 620. / 27.;
    cache_ps[index_p(cc, 3, 4, 1)] = 268. / 27. * L_1 - 64. / 9. * L_2 - 16. / 3. * z - 16. / 9. * M_PI2 + 112. / (9. * sqrt3) * M_PI + 3506. / 81.;
    cache_ps[index_p(cc, 3, 5, 1)] = 5632. / 9. * L_1 - 4096. / 3. * L_2 - 8640. * z - 256. / 3. * M_PI2 + 256. / (3. * sqrt3) * M_PI + 9728. / 27.;
    cache_ps[index_p(cc, 3, 6, 1)] = 9184. / 27. * L_1 - 1024. / 9. * L_2 - 160. / 3. * z - 256. / 9. * M_PI2 + 1072. / (9. * sqrt3) * M_PI + 88688. / 81.;
    cache_ps[index_p(cc, 4, 4, 1)] = 1028. / 81. * L_1 + 136. / 27. * L_2 - 8. / 3. * M_PI2 * z + 230. / 9. * z - 134. / 27. * M_PI2 + 230. / (27. * sqrt3) * M_PI + 6214. / 243.;
    cache_ps[index_p(cc, 4, 5, 1)] = 9472. / 27. * L_1 - 1024. / 9. * L_2 + 608. / 3. * z - 256. / 9. * M_PI2 + 1792. / (9. * sqrt3) + 64784. / 81.;
    cache_ps[index_p(cc, 4, 6, 1)] = 10792. / 81. * L_1 + 2048. / 27. * L_2 - 160. / 3. * M_PI2 * z + 3568. / 9. * z - 2560. / 27. * M_PI2 + 4264. / (27. * sqrt3) * M_PI + 123080. / 243.;
    cache_ps[index_p(cc, 5, 5, 1)] = 45056. / 9. * L_1 - 28160. / 3. * L_2 - 58752. * z - 2048. / 3. * M_PI2 - 2048. / (3. * sqrt3) * M_PI - 349184. / 27.;
    cache_ps[index_p(cc, 5, 6, 1)] = 167680. / 27. * L_1 - 16384. / 9. * L_2 + 6080. / 3. * z - 4096. / 9. * M_PI2 + 17152. / (9. * sqrt3) * M_PI + 1502720. / 81.;
    cache_ps[index_p(cc, 6, 6, 1)] = 75392. / 81. * L_1 + 11776. / 27. * L_2 - 1088. / 3. * M_PI2 * z + 23696. / 9. * z - 17024. / 27. * M_PI2 + 19712. / (27. * sqrt3) * M_PI + 717184. / 243.;
 
    // Corrections from LO coefficients when resummation logz
    cache_p[index_p(cc, 5, 5, 1)] += -1296. * n_v * z_replace;
    cache_p[index_p(cc, 6, 6, 1)] += -72. * n_v * z_replace;
 
    // Gerlach thesis eq. (6.17)
    for (int i = 3; i <= 6; i++)
    {
        for (int j = i; j <= 6; j++)
        {
            cache_p[index_p(cu, i, j, 1)] = cache_p[index_p(cc, i, j, 1)];
            cache_p[index_p(uu, i, j, 1)] = cache_p[index_p(cc, i, j, 1)];
            cache_ps[index_p(cu, i, j, 1)] = cache_ps[index_p(cc, i, j, 1)];
            cache_ps[index_p(uu, i, j, 1)] = cache_ps[index_p(cc, i, j, 1)];
        }
    }
 
    // Gerlach thesis eq. (6.19)
    cache_p[index_p(cc, 1, 8, 1)] = 5. / 18. * sqrt1minus4z + 5. / 9. * z * sqrt1minus4z;
    cache_p[index_p(cc, 2, 8, 1)] = -5. / 3. * sqrt1minus4z - 10. / 3. * z * sqrt1minus4z;
    cache_ps[index_p(cc, 1, 8, 1)] = 4. / 9. * sqrt1minus4z + 8. / 9. * z * sqrt1minus4z;
    cache_ps[index_p(cc, 2, 8, 1)] = -8. / 3. * sqrt1minus4z - 16. / 3. * z * sqrt1minus4z;
    cache_p[index_p(uu, 1, 8, 1)] = 5. / 18.;
    cache_p[index_p(uu, 2, 8, 1)] = -5. / 3.;
    cache_ps[index_p(uu, 1, 8, 1)] = 4. / 9.;
    cache_ps[index_p(uu, 2, 8, 1)] = -8. / 3.;
 
    // Gerlach thesis eq. (6.27)
    for (int i = 1; i <= 2; i++)
    {
        cache_p[index_p(cu, i, 8, 1)] = 0.5 * (cache_p[index_p(cc, i, 8, 1)] + cache_p[index_p(uu, i, 8, 1)]);
        cache_ps[index_p(cu, i, 8, 1)] = 0.5 * (cache_ps[index_p(cc, i, 8, 1)] + cache_ps[index_p(uu, i, 8, 1)]);
    }
 
    // Gerlach thesis eq. (6.21)
    cache_p[index_p(cc, 3, 8, 1)] = -32. / 3.;
    cache_p[index_p(cc, 4, 8, 1)] = (-139. - 30. * sqrt1minus4z - 60. * sqrt1minus4z * z) / 18.;
    cache_p[index_p(cc, 5, 8, 1)] = -512. / 3.;
    cache_p[index_p(cc, 6, 8, 1)] = (-1684. - 300. * sqrt1minus4z - 600. * sqrt1minus4z * z) / 18.;
    cache_ps[index_p(cc, 3, 8, 1)] = 64. / 3.;
    cache_ps[index_p(cc, 4, 8, 1)] = (4. * (1. - 6. * sqrt1minus4z - 12. * sqrt1minus4z * z)) / 9.;
    cache_ps[index_p(cc, 5, 8, 1)] = 1024. / 3.;
    cache_ps[index_p(cc, 6, 8, 1)] = (4 * (124. - 60. * sqrt1minus4z - 120. * sqrt1minus4z * z)) / 9.;
 
    // Gerlach thesis eq. (6.20)
    for (int i = 3; i <= 6; i++)
    {
        cache_p[index_p(cu, i, 8, 1)] = cache_p[index_p(cc, i, 8, 1)];
        cache_p[index_p(uu, i, 8, 1)] = cache_p[index_p(cc, i, 8, 1)];
        cache_ps[index_p(cu, i, 8, 1)] = cache_ps[index_p(cc, i, 8, 1)];
        cache_ps[index_p(uu, i, 8, 1)] = cache_ps[index_p(cc, i, 8, 1)];
    }
 
    double L_12 = L_1 * L_1;
    double L_22 = L_2 * L_2;
 
    gslpp::complex log1M1OverSqrtz = log(abs(1. - 1. / sqrtz));
    if (1. / sqrtz > 1)
        log1M1OverSqrtz += M_PI * complex::i();
    gslpp::complex logM1OverSqrtz = log(abs(-1. / sqrtz)) + M_PI * complex::i();
    gslpp::complex logM1OverSqrtz2 = logM1OverSqrtz * logM1OverSqrtz;
 
    // LO in z
    cache_p[index_p(cc, 1, 1, 2)] = (814589597. / 4199040. + (1320817. * L_1) / 17496. + (3211. * L_12) / 324. + (259603. * L_2) / 11664. + (12911. * L_1 * L_2) / 972. - (311. * L_22) / 216. + (56. * Cl2PI3 * sqrt3) / 729. -
                                     ((5. * complex::i()) / 3888.) * log3 * log3 * sqrt3 + (855371. * sqrtz) / 180000. + (22. * log2z * sqrtz) / 9. + (88. * logM1OverSqrtz2 * sqrtz) / 9. + (88. * logM1OverSqrtz * logz * sqrtz) / 9. -
                                     ((5. * -1.) / 486.) * sqrt3 * t_2 - (22. * log2z) / (9. * z) - (88. * logM1OverSqrtz2) / (9. * z) - (88. * logM1OverSqrtz * logz) / (9. * z) + (22. * log2z) / (9. * sqrtz) +
                                     (88. * logM1OverSqrtz2) / (9. * sqrtz) + (88. * logM1OverSqrtz * logz) / (9. * sqrtz) - (3547371902315179. * z) / 3455518320000. - (721919. * L_1 * z) / 1944. - (2347. * L_12 * z) / 54. +
                                     (1891. * L_2 * z) / 81. - (337. * L_1 * L_2 * z) / 9. + (187. * L_22 * z) / 18. - (88. * logM1OverSqrtz2 * z) / 9. - (4823. * logz * z) / 9. - (1348. * L_1 * logz * z) / 9. - (88. * L_2 * logz * z) / 3. -
                                     (88. * logM1OverSqrtz * logz * z) / 9. - (28333. * zeta3) / 486. + (128581. * z * zeta3) / 216. - (5. * sqrt3 * (20. + 30. * L_1 + 3. * log3 + 180. * z) * (M_PI)) / 17496. +
                                     ((-684288. + 684288. * sqrtz - (216641. + 87480. * L_1 - 432. * L_2 + 146016. * log12sqrt52 + 5112. * log2 - (50. * complex::i()) * sqrt3 + 158184. * sqrt5 - 684288. * sqrtz) * z +
                                       18. * (197453. + 2232. * L_1 + 912. * L_2 - 9792. * log12sqrt52 - 26508. * log2 + 576. * logz + 9744. * sqrt5) * z2) *
                                      (M_PI2)) /
                                         (69984. * z) +
                                     (23. * (M_PI4)) / 4860. - (13637. * z * (M_PI4)) / 116640.)
                                        .real();
    cache_p[index_p(cc, 1, 2, 2)] = (-95740679. / 349920. - (1026907. * L_1) / 5832. - (1751. * L_12) / 54. - (619. * L_2) / 972. + (166. * L_1 * L_2) / 81. - L_22 / 18. - (224. * Cl2PI3 * sqrt3) / 243. + ((5. * complex::i()) / 324.) * log3 * log3 * sqrt3 +
                                     (2895091. * sqrtz) / 60000. + (8. * log2z * sqrtz) / 3. + (32. * logM1OverSqrtz2 * sqrtz) / 3. + (32. * logM1OverSqrtz * logz * sqrtz) / 3. + ((10. * -1.) / 81.) * sqrt3 * t_2 - (8. * log2z) / (3. * z) -
                                     (32. * logM1OverSqrtz2) / (3. * z) - (32. * logM1OverSqrtz * logz) / (3. * z) + (8. * log2z) / (3. * sqrtz) + (32. * logM1OverSqrtz2) / (3. * sqrtz) +
                                     (32. * logM1OverSqrtz * logz) / (3. * sqrtz) + (2370097879. * z) / 1944000. + (117443. * L_1 * z) / 162. + (1193. * L_12 * z) / 9. + (6350. * L_2 * z) / 27. + (64. * L_1 * L_2 * z) / 3. + (34. * L_22 * z) / 3. -
                                     (32. * logM1OverSqrtz2 * z) / 3. + 124. * logz * z + (256. * L_1 * logz * z) / 3. - 32. * L_2 * logz * z - (32. * logM1OverSqrtz * logz * z) / 3. - (10573. * zeta3) / 324. + (85027. * z * zeta3) / 90. +
                                     (5. * sqrt3 * (20. + 30. * L_1 + 3. * log3 + 180. * z) * (M_PI)) / 1458. -
                                     ((622080. - 622080. * sqrtz - 5. * (497221. + 116640. * L_1 - 864. * L_2 - 39744. * log12sqrt52 + 85824. * log2 - (100. * complex::i()) * sqrt3 - 43056. * sqrt5 + 124416. * sqrtz) * z +
                                       72. * (64865. + 3960. * L_1 + 2280. * L_2 + 3168. * log12sqrt52 + 35070. * log2 + 1440. * logz - 3288. * sqrt5) * z2) *
                                      (M_PI2)) /
                                         (58320. * z) -
                                     (799. * (M_PI4)) / 810. + (20833. * z * (M_PI4)) / 4860.)
                                        .real();
    cache_p[index_p(cc, 2, 2, 2)] = (74041. / 14580. + (106199. * L_1) / 972. + (239. * L_12) / 18. - (5117. * L_2) / 81. - (202. * L_1 * L_2) / 27. - (19. * L_22) / 3. + (224. * Cl2PI3 * sqrt3) / 81. - ((5. * complex::i()) / 108.) * log3 * log3 * sqrt3 -
                                     (1052759. * sqrtz) / 10000. + 4. * log2z * sqrtz + 16. * logM1OverSqrtz2 * sqrtz + 16. * logM1OverSqrtz * logz * sqrtz - ((10. * -1.) / 27.) * sqrt3 * t_2 - (4. * log2z) / z -
                                     (16. * logM1OverSqrtz2) / z - (16. * logM1OverSqrtz * logz) / z + (4. * log2z) / sqrtz + (16. * logM1OverSqrtz2) / sqrtz + (16. * logM1OverSqrtz * logz) / sqrtz -
                                     (69930883. * z) / 101250. - (16463. * L_1 * z) / 54. - (122. * L_12 * z) / 3. - (109. * L_2 * z) / 18. - 22. * L_1 * L_2 * z + 17. * L_22 * z - 16. * logM1OverSqrtz2 * z - 496. * logz * z - 88. * L_1 * logz * z -
                                     48. * L_2 * logz * z - 16. * logM1OverSqrtz * logz * z - (3157. * zeta3) / 54. + (2521. * z * zeta3) / 15. - (5. * sqrt3 * (20. + 30. * L_1 + 3. * log3 + 180. * z) * (M_PI)) / 486. +
                                     (-177247. / 3888. - (4. * log12sqrt52) / 9. + (148. * log2) / 27. + ((25. * complex::i()) / 972.) * sqrt3 - (13. * sqrt5) / 27. + 16. * sqrtz - 16. / z + 16. / sqrtz + (5369. * z) / 108. - (16. * log12sqrt52 * z) / 15. -
                                      (178. * log2 * z) / 9. + (16. * logz * z) / 3. + (28. * sqrt5 * z) / 45. + L_1 * (-15. + (26. * z) / 3.) + (2. * L_2 * (1. + 38. * z)) / 9.) *
                                         (M_PI2) +
                                     (971. * (M_PI4)) / 540. + (5203. * z * (M_PI4)) / 3240.)
                                        .real();
 
    cache_ps[index_p(cc, 1, 1, 2)] = (-67489177. / 262440. - (77617. * L_1) / 2187. - (902. * L_12) / 243. -
                                      (5504. * L_2) / 729. - (3064. * L_1 * L_2) / 243. + (260. * L_22) / 27. +
                                      (104. * Cl2PI3 * sqrt3) / 729. - (complex::i() / 486.) * log3 * log3 * sqrt3 +
                                      (166687. * sqrtz) / 6000. - ((4. * -1.) / 243.) * sqrt3 * t_2 -
                                      (27341897029. * z) / 29160000. - (97999. * L_1 * z) / 243. - (9272. * L_2 * z) / 81. -
                                      (9272. * logz * z) / 27. + (28528. * zeta3) / 243. + (29. * z * zeta3) / 3. -
                                      (sqrt3 * (20. + 30. * L_1 + 3. * log3 + 180. * z) * (M_PI)) / 2187. -
                                      ((44209. - 37152. * log12sqrt52 + 126216. * log2 - (10. * complex::i()) * sqrt3 -
                                        40248. * sqrt5 + 127854. * z - 197208. * log2 * z + 10368. * logz * z +
                                        37152. * sqrt5 * z + 17496. * L_1 * (1. + 2. * z) + 1728. * L_2 * (1. + 2. * z)) *
                                       (M_PI2)) /
                                          8748. +
                                      (449. * (M_PI4)) / 1215. - (27529. * z * (M_PI4)) / 14580.)
                                         .real();
    cache_ps[index_p(cc, 1, 2, 2)] = (-1336127. / 2187. - (28733. * L_1) / 729. + (44. * L_12) / 81. - (2176. * L_2) / 243. -
                                      (1856. * L_1 * L_2) / 81. + (208. * L_22) / 9. - (416. * Cl2PI3 * sqrt3) / 243. +
                                      ((2. * complex::i()) / 81.) * log3 * log3 * sqrt3 + (8773. * sqrtz) / 100. + ((16. * -1.) / 81.) * sqrt3 * t_2 +
                                      (1672496939. * z) / 4860000. - (23372. * L_1 * z) / 81. - (5248. * L_2 * z) / 27. -
                                      (5248. * logz * z) / 9. + (13934. * zeta3) / 81. - 244. * z * zeta3 +
                                      (4. * sqrt3 * (20. + 30. * L_1 + 3. * log3 + 180. * z) * (M_PI)) / 729. +
                                      ((39995. + 8640. * log12sqrt52 - 29232. * log2 - (20. * complex::i()) * sqrt3 + 9360. * sqrt5 +
                                        58284. * z + 11232. * log2 * z + 20736. * logz * z - 8640. * sqrt5 * z +
                                        23328. * L_1 * (1. + 2. * z) + 3456. * L_2 * (1. + 2. * z)) *
                                       (M_PI2)) /
                                          1458. +
                                      (226. * (M_PI4)) / 405. + (5692. * z * (M_PI4)) / 1215.)
                                         .real();
    cache_ps[index_p(cc, 2, 2, 2)] = (27476329. / 58320. + (40370. * L_1) / 243. + (604. * L_12) / 27. + (6928. * L_2) / 81. +
                                      (1064. * L_1 * L_2) / 27. - (52. * L_22) / 3. + (416. * Cl2PI3 * sqrt3) / 81. -
                                      ((2. * complex::i()) / 27.) * log3 * log3 * sqrt3 - (26319. * sqrtz) / 250. - ((16. * -1.) / 27.) * sqrt3 * t_2 +
                                      (287805209. * z) / 81000. + (18836. * L_1 * z) / 27. + (928. * L_2 * z) / 9. + (928. * logz * z) / 3. -
                                      (4388. * zeta3) / 27. - 600. * z * zeta3 -
                                      (4. * sqrt3 * (20. + 30. * L_1 + 3. * log3 + 180. * z) * (M_PI)) / 243. -
                                      ((41279. - 1728. * log12sqrt52 + 11808. * log2 - (20. * complex::i()) * sqrt3 - 1872. * sqrt5 +
                                        144684. * z - 14688. * log2 * z + 20736. * logz * z + 1728. * sqrt5 * z +
                                        11664. * L_1 * (1. + 2. * z) + 3456. * L_2 * (1. + 2. * z)) *
                                       (M_PI2)) /
                                          486. +
                                      (398. * (M_PI4)) / 135. + (7991. * z * (M_PI4)) / 405.)
                                         .real();
 
    cache_p[index_p(uu, 1, 1, 2)] = (814589597. / 4199040. + (1320817. * L_1) / 17496. + (3211. * L_12) / 324. +
                                     (259603. * L_2) / 11664. + (12911. * L_1 * L_2) / 972. - (311. * L_22) / 216. +
                                     (56. * Cl2PI3 * sqrt3) / 729. - ((5. * complex::i()) / 3888.) * log3 * log3 * sqrt3 +
                                     (855371. * sqrtz) / 180000. - ((5. * -1.) / 486.) * sqrt3 * t_2 +
                                     (5298817581707. * z) / 1151839440000. + (5. * L_1 * z) / 81. - (50. * logz * z) / 9. -
                                     (28333. * zeta3) / 486. - (5. * (20. + 30. * L_1 + 3. * log3) * sqrt3 * (M_PI)) / 17496. -
                                     ((216641. + 87480. * L_1 - 432. * L_2 + 146016. * log12sqrt52 + 5112. * log2 -
                                       (50. * complex::i()) * sqrt3 + 158184. * sqrt5) *
                                      (M_PI2)) /
                                         69984. +
                                     (23. * (M_PI4)) / 4860.)
                                        .real();
    cache_p[index_p(uu, 1, 2, 2)] = (-95740679. / 349920. - (1026907. * L_1) / 5832. - (1751. * L_12) / 54. - (619. * L_2) / 972. +
                                     (166. * L_1 * L_2) / 81. - L_22 / 18. - (224. * Cl2PI3 * sqrt3) / 243. +
                                     ((5. * complex::i()) / 324.) * log3 * log3 * sqrt3 + (2895091. * sqrtz) / 60000. +
                                     ((10. * -1.) / 81.) * sqrt3 * t_2 - (55644107. * z) / 648000. - (20. * L_1 * z) / 27. +
                                     (8. * logz * z) / 3. - (10573. * zeta3) / 324. +
                                     (5. * (20. + 30. * L_1 + 3. * log3) * sqrt3 * (M_PI)) / 1458. +
                                     ((497221. + 116640. * L_1 - 864. * L_2 - 39744. * log12sqrt52 + 85824. * log2 -
                                       (100. * complex::i()) * sqrt3 - 43056. * sqrt5) *
                                      (M_PI2)) /
                                         11664. -
                                     (799. * (M_PI4)) / 810.)
                                        .real();
    cache_p[index_p(uu, 2, 2, 2)] = (74041. / 14580. + (106199. * L_1) / 972. + (239. * L_12) / 18. - (5117. * L_2) / 81. -
                                     (202. * L_1 * L_2) / 27. - (19. * L_22) / 3. + (224. * Cl2PI3 * sqrt3) / 81. -
                                     ((5. * complex::i()) / 108.) * log3 * log3 * sqrt3 - (1052759. * sqrtz) / 10000. -
                                     ((10. * -1.) / 27.) * sqrt3 * t_2 + (9532631. * z) / 135000. + (20. * L_1 * z) / 9. - 32. * logz * z -
                                     (3157. * zeta3) / 54. - (5. * (20. + 30. * L_1 + 3. * log3) * sqrt3 * (M_PI)) / 486. -
                                     ((177247. + 58320. * L_1 - 864. * L_2 + 1728. * log12sqrt52 - 21312. * log2 -
                                       (100. * complex::i()) * sqrt3 + 1872. * sqrt5) *
                                      (M_PI2)) /
                                         3888. +
                                     (971. * (M_PI4)) / 540.)
                                        .real();
    cache_ps[index_p(uu, 1, 1, 2)] = (-67489177. / 262440. - (77617. * L_1) / 2187. - (902. * L_12) / 243. - (5504. * L_2) / 729. -
                                      (3064. * L_1 * L_2) / 243. + (260. * L_22) / 27. + (104. * Cl2PI3 * sqrt3) / 729. -
                                      (complex::i() / 486.) * log3 * log3 * sqrt3 + (166687. * sqrtz) / 6000. - ((4. * -1.) / 243.) * sqrt3 * t_2 -
                                      (261095543. * z) / 9720000. + (8. * L_1 * z) / 81. + (28528. * zeta3) / 243. -
                                      ((20. + 30. * L_1 + 3. * log3) * sqrt3 * (M_PI)) / 2187. -
                                      ((44209. + 17496. * L_1 + 1728. * L_2 - 37152. * log12sqrt52 + 126216. * log2 -
                                        (10. * complex::i()) * sqrt3 - 40248. * sqrt5) *
                                       (M_PI2)) /
                                          8748. +
                                      (449. * (M_PI4)) / 1215.)
                                         .real();
    cache_ps[index_p(uu, 1, 2, 2)] = (-1336127. / 2187. - (28733. * L_1) / 729. + (44. * L_12) / 81. - (2176. * L_2) / 243. -
                                      (1856. * L_1 * L_2) / 81. + (208. * L_22) / 9. - (416. * Cl2PI3 * sqrt3) / 243. +
                                      ((2. * complex::i()) / 81.) * log3 * log3 * sqrt3 + (8773. * sqrtz) / 100. + ((16. * -1.) / 81.) * sqrt3 * t_2 -
                                      (117168887. * z) / 1620000. - (32. * L_1 * z) / 27. + (13934. * zeta3) / 81. +
                                      (4. * (20. + 30. * L_1 + 3. * log3) * sqrt3 * (M_PI)) / 729. +
                                      ((39995. + 23328. * L_1 + 3456. * L_2 + 8640. * log12sqrt52 - 29232. * log2 -
                                        (20. * complex::i()) * sqrt3 + 9360. * sqrt5) *
                                       (M_PI2)) /
                                          1458. +
                                      (226. * (M_PI4)) / 405.)
                                         .real();
    cache_ps[index_p(uu, 2, 2, 2)] = (27476329. / 58320. + (40370. * L_1) / 243. + (604. * L_12) / 27. + (6928. * L_2) / 81. +
                                      (1064. * L_1 * L_2) / 27. - (52. * L_22) / 3. + (416. * Cl2PI3 * sqrt3) / 81. -
                                      ((2. * complex::i()) / 27.) * log3 * log3 * sqrt3 - (26319. * sqrtz) / 250. - ((16. * -1.) / 27.) * sqrt3 * t_2 +
                                      (4778443. * z) / 27000. + (32. * L_1 * z) / 9. - (4388. * zeta3) / 27. -
                                      (4. * (20. + 30. * L_1 + 3. * log3) * sqrt3 * (M_PI)) / 243. -
                                      ((41279. + 11664. * L_1 + 3456. * L_2 - 1728. * log12sqrt52 + 11808. * log2 -
                                        (20. * complex::i()) * sqrt3 - 1872. * sqrt5) *
                                       (M_PI2)) /
                                          486. +
                                      (398. * (M_PI4)) / 135.)
                                         .real();
 
    // Corrections from LO and NLO coefficients (incl. pengFlag)
    // resummation of logz
    cache_p[index_p(cc, 1, 1, 2)] += (-2. * (1661. + 1113. * log2z + 132. * M_PI2) * z + logz * (-30. + (39727. + 12132. * L_1 + 2376. * L_2 - 12. * M_PI2) * z)) / 81.;
    cache_p[index_p(cc, 1, 2, 2)] += (-8. * ((151 + 147. * log2z + 12. * M_PI2) * z + 2 * logz * (-12. + (448. + 144. * L_1 - 54. * L_2 - 3. * M_PI2) * z))) / 27.;
    cache_p[index_p(cc, 2, 2, 2)] += (-4. * (66. * logz - 2. * logz * (518. + 99. * L_1 + 54. * L_2 - 6. * M_PI2) * z + (151. + 21. * log2z + 12. * M_PI2) * z)) / 9.;
    cache_ps[index_p(cc, 1, 1, 2)] += (8. * logz * (4. + (507. + 172. * logz + 4. * M_PI2) * z)) / 27.;
    cache_ps[index_p(cc, 1, 2, 2)] += (32. * logz * (-1. + 4. * (12. + 5. * logz - M_PI2) * z)) / 9.;
    cache_ps[index_p(cc, 2, 2, 2)] += (32. * logz * (-2. + (-9. + 4. * logz + 4. * M_PI2) * z)) / 3.;
 
    // resummation affecting arguments of logs and Dilogs
    cache_p[index_p(cc, 1, 1, 2)] += 2. / 81. * logz * (15. + 1615. * z + 1014. * z * logz);
    cache_p[index_p(cc, 1, 2, 2)] += 4. / 27. * logz * (-48. + 973. * z + 276. * z * logz);
    cache_p[index_p(cc, 2, 2, 2)] += 8. / 9. * logz * (33. + 4. * z + 6. * z * logz);
    cache_ps[index_p(cc, 1, 1, 2)] += -32. / 27. * logz + 5216. / 27. * logz * z - 1376. / 27. * log2z * z;
    cache_ps[index_p(cc, 1, 2, 2)] += 32. / 9. * logz + 3712. / 9. * logz * z - 640. / 9. * log2z * z;
    cache_ps[index_p(cc, 2, 2, 2)] += 64. / 3. * logz - 640. / 3. * logz * z - 128. / 3. * log2z * z;
 
    double L_b = 2. * log(mu_b / Mb);
 
    // Transforming mb in NLO coefficients from Pole scheme to MSbar
    cache_p[index_p(cc, 1, 1, 2)] += (8. * (4. + 3. * L_b) * (-196. + 675. * z)) / 243.;
    cache_p[index_p(cc, 1, 2, 2)] += (-44. * (4. + 3. * L_b) * (-19. + 108. * z)) / 81.;
    cache_p[index_p(cc, 2, 2, 2)] += (-16. * (4. + 3. * L_b) * (-4. + 27. * z)) / 27.;
    cache_ps[index_p(cc, 1, 1, 2)] += (1264. * (4. + 3. * L_b)) / 243.;
    cache_ps[index_p(cc, 1, 2, 2)] += (416. * (4. + 3. * L_b)) / 81.;
    cache_ps[index_p(cc, 2, 2, 2)] += (-704. * (4. + 3. * L_b)) / 27.;
 
    cache_p[index_p(uu, 1, 1, 2)] += (-1568. * (4. + 3. * L_1)) / 243.;
    cache_p[index_p(uu, 1, 2, 2)] += (836. * (4. + 3. * L_1)) / 81.;
    cache_p[index_p(uu, 2, 2, 2)] += (64. * (4. + 3. * L_1)) / 27.;
    cache_ps[index_p(uu, 1, 1, 2)] += (1264. * (4. + 3. * L_b)) / 243.;
    cache_ps[index_p(uu, 1, 2, 2)] += (416. * (4. + 3. * L_b)) / 81.;
    cache_ps[index_p(uu, 2, 2, 2)] += (-704. * (4. + 3. * L_b)) / 27.;
 
    // Gerlach thesis eq. (6.27)
    for (int i = 1; i <= 2; i++)
    {
        for (int j = i; j <= 2; j++)
        {
            cache_p[index_p(cu, i, j, 2)] = 0.5 * (cache_p[index_p(cc, i, j, 2)] + cache_p[index_p(uu, i, j, 2)]);
            cache_ps[index_p(cu, i, j, 2)] = 0.5 * (cache_ps[index_p(cc, i, j, 2)] + cache_ps[index_p(uu, i, j, 2)]);
        }
    }
 
    // LO in z
    cache_p[index_p(cu, 1, 1, 2)] = (814589597. / 4199040. + (56. * Cl2PI3) / (243. * sqrt3) + (1320817. * L_1) / 17496. + (3211. * L_12) / 324. +
                                     (259603. * L_2) / 11664. + (12911. * L_1 * L_2) / 972. - (311. * L_22) / 216. - (((5. * complex::i()) / 1296.) * log3 * log3) / sqrt3 +
                                     (855371. * sqrtz) / 180000. + (11. * log2z * sqrtz) / 9. + (44. * logM1OverSqrtz2 * sqrtz) / 9. +
                                     (44. * logM1OverSqrtz * logz * sqrtz) / 9. - (((5. * -1.) / 162.) * t_2) / sqrt3 - (11. * log2z) / (9. * z) -
                                     (44. * logM1OverSqrtz2) / (9. * z) - (44. * logM1OverSqrtz * logz) / (9. * z) + (11. * log2z) / (9. * sqrtz) +
                                     (44. * logM1OverSqrtz2) / (9. * sqrtz) + (44. * logM1OverSqrtz * logz) / (9. * sqrtz) -
                                     (1765737724785029. * z) / 3455518320000. - (721799. * L_1 * z) / 3888. - (2347. * L_12 * z) / 108. + (1891. * L_2 * z) / 162. -
                                     (337. * L_1 * L_2 * z) / 18. + (187. * L_22 * z) / 36. - (44. * logM1OverSqrtz2 * z) / 9. - (4873. * logz * z) / 18. -
                                     (674. * L_1 * logz * z) / 9. - (44. * L_2 * logz * z) / 3. - (44. * logM1OverSqrtz * logz * z) / 9. - (28333. * zeta3) / 486. +
                                     (128581. * z * zeta3) / 432. - (5. * (3. * sqrt3 * log3 + 30. * L_1 * sqrt3 + 10. * sqrt3 * (2. + 9. * z)) * (M_PI)) / 17496. +
                                     ((-342144. + 342144. * sqrtz + (-216641. + (50. * complex::i()) * sqrt3 - 87480. * L_1 + 432. * L_2 - 5112. * log2 - 158184. * sqrt5 + 342144. * sqrtz) * z + 9. * (197453. + 2232. * L_1 + 912. * L_2 - 26508. * log2 + 576. * logz + 9744. * sqrt5) * z2 - 864. * z * (169. + 102. * z) * log12sqrt52) * (M_PI2)) / (69984. * z) +
                                     (23. * (M_PI4)) / 4860. - (13637. * z * (M_PI4)) / 233280.)
                                        .real();
    cache_p[index_p(cu, 1, 2, 2)] = (-95740679. / 349920. - (1026907. * L_1) / 5832. - (1751. * L_12) / 54. - (619. * L_2) / 972. + (166. * L_1 * L_2) / 81. - L_22 / 18. -
                                     (224. * Cl2PI3 * sqrt3) / 243. + ((5. * complex::i()) / 324.) * log3 * log3 * sqrt3 + (2895091. * sqrtz) / 60000. + (4. * log2z * sqrtz) / 3. +
                                     (16. * logM1OverSqrtz2 * sqrtz) / 3. + (16. * logM1OverSqrtz * logz * sqrtz) / 3. + ((10. * -1.) / 81.) * sqrt3 * t_2 -
                                     (4. * log2z) / (3. * z) - (16. * logM1OverSqrtz2) / (3. * z) - (16. * logM1OverSqrtz * logz) / (3. * z) +
                                     (4. * log2z) / (3. * sqrtz) + (16. * logM1OverSqrtz2) / (3. * sqrtz) +
                                     (16. * logM1OverSqrtz * logz) / (3. * sqrtz) + (1101582779. * z) / 1944000. + (117323. * L_1 * z) / 324. +
                                     (1193. * L_12 * z) / 18. + (3175. * L_2 * z) / 27. + (32. * L_1 * L_2 * z) / 3. + (17. * L_22 * z) / 3. - (16. * logM1OverSqrtz2 * z) / 3. +
                                     (190. * logz * z) / 3. + (128. * L_1 * logz * z) / 3. - 16. * L_2 * logz * z - (16. * logM1OverSqrtz * logz * z) / 3. -
                                     (10573. * zeta3) / 324. + (85027. * z * zeta3) / 180. + (5. * sqrt3 * (20. + 30. * L_1 + 3. * log3 + 90. * z) * (M_PI)) / 1458. -
                                     ((311040. - 311040. * sqrtz - 5. * (497221. + 116640. * L_1 - 864. * L_2 - 39744. * log12sqrt52 + 85824. * log2 - (100. * complex::i()) * sqrt3 - 43056. * sqrt5 + 62208. * sqrtz) * z +
                                       36. * (64865. + 3960. * L_1 + 2280. * L_2 + 3168. * log12sqrt52 + 35070. * log2 + 1440. * logz - 3288. * sqrt5) * z2) *
                                      (M_PI2)) /
                                         (58320. * z) -
                                     (799. * (M_PI4)) / 810. + (20833. * z * (M_PI4)) / 9720.)
                                        .real();
    cache_p[index_p(cu, 2, 2, 2)] = (74041. / 14580. + (106199. * L_1) / 972. + (239. * L_12) / 18. - (5117. * L_2) / 81. - (202. * L_1 * L_2) / 27. - (19. * L_22) / 3. +
                                     (224. * Cl2PI3 * sqrt3) / 81. - ((5. * complex::i()) / 108.) * log3 * log3 * sqrt3 - (1052759. * sqrtz) / 10000. + 2. * log2z * sqrtz +
                                     8. * logM1OverSqrtz2 * sqrtz + 8. * logM1OverSqrtz * logz * sqrtz - ((10. * -1.) / 27.) * sqrt3 * t_2 - (2. * log2z) / z -
                                     (8. * logM1OverSqrtz2) / z - (8. * logM1OverSqrtz * logz) / z + (2. * log2z) / sqrtz +
                                     (8. * logM1OverSqrtz2) / sqrtz + (8. * logM1OverSqrtz * logz) / sqrtz - (251125639. * z) / 810000. -
                                     (16343. * L_1 * z) / 108. - (61. * L_12 * z) / 3. - (109. * L_2 * z) / 36. - 11. * L_1 * L_2 * z + (17. * L_22 * z) / 2. -
                                     8. * logM1OverSqrtz2 * z - 264. * logz * z - 44. * L_1 * logz * z - 24. * L_2 * logz * z - 8. * logM1OverSqrtz * logz * z -
                                     (3157. * zeta3) / 54. + (2521. * z * zeta3) / 30. - (5. * sqrt3 * (20. + 30. * L_1 + 3. * log3 + 90. * z) * (M_PI)) / 486. +
                                     (-177247. / 3888. - (4. * log12sqrt52) / 9. + (148. * log2) / 27. + ((25. * complex::i()) / 972.) * sqrt3 - (13. * sqrt5) / 27. +
                                      8. * sqrtz - 8. / z + 8. / sqrtz + (5369. * z) / 216. - (8. * log12sqrt52 * z) / 15. - (89. * log2 * z) / 9. +
                                      (8. * logz * z) / 3. + (14. * sqrt5 * z) / 45. + L_1 * (-15. + (13. * z) / 3.) + (2. * L_2 * (1. + 19. * z)) / 9.) *
                                         (M_PI2) +
                                     (971. * (M_PI4)) / 540. + (5203. * z * (M_PI4)) / 6480.)
                                        .real();
    cache_ps[index_p(cu, 1, 1, 2)] = (-67489177. / 262440. + (104. * Cl2PI3) / (243. * sqrt3) - (77617. * L_1) / 2187. -
                                      (902. * L_12) / 243. - (5504. * L_2) / 729. - (3064. * L_1 * L_2) / 243. + (260. * L_22) / 27. -
                                      ((complex::i() / 162.) * log3 * log3) / sqrt3 + (166687. * sqrtz) / 6000. - (((4. * -1.) / 81.) * t_2) / sqrt3 -
                                      (14062591829. * z) / 29160000. - (97975. * L_1 * z) / 486. - (4636. * L_2 * z) / 81. -
                                      (4636. * logz * z) / 27. + (28528. * zeta3) / 243. + (29. * z * zeta3) / 6. -
                                      ((3. * sqrt3 * log3 + 30. * L_1 * sqrt3 + 10. * sqrt3 * (2. + 9. * z)) * (M_PI)) / 2187. +
                                      ((-44209. + (10. * complex::i()) * sqrt3 - 1728. * L_2 - 163368. * log2 + 40248. * sqrt5 - 63927. * z -
                                        1728. * L_2 * z + 98604. * log2 * z - 5184. * logz * z - 18576. * sqrt5 * z - 17496. * L_1 * (1. + z) +
                                        37152. * (log12sqrt52 + log2)) *
                                       (M_PI2)) /
                                          8748. +
                                      (449. * (M_PI4)) / 1215. -
                                      (27529. * z * (M_PI4)) / 29160.)
                                         .real();
    cache_ps[index_p(cu, 1, 2, 2)] = (-1336127. / 2187. - (28733. * L_1) / 729. + (44. * L_12) / 81. - (2176. * L_2) / 243. - (1856. * L_1 * L_2) / 81. +
                                      (208. * L_22) / 9. - (416. * Cl2PI3 * sqrt3) / 243. + ((2. * complex::i()) / 81.) * log3 * log3 * sqrt3 +
                                      (8773. * sqrtz) / 100. + ((16. * -1.) / 81.) * sqrt3 * t_2 + (660495139. * z) / 4860000. -
                                      (11734. * L_1 * z) / 81. - (2624. * L_2 * z) / 27. - (2624. * logz * z) / 9. + (13934. * zeta3) / 81. -
                                      122. * z * zeta3 + (4. * sqrt3 * (20. + 30. * L_1 + 3. * log3 + 90. * z) * (M_PI)) / 729. +
                                      ((39995. + 8640. * log12sqrt52 - 29232. * log2 - (20. * complex::i()) * sqrt3 + 9360. * sqrt5 + 29142. * z +
                                        5616. * log2 * z + 10368. * logz * z - 4320. * sqrt5 * z + 23328. * L_1 * (1. + z) +
                                        3456. * L_2 * (1. + z)) *
                                       (M_PI2)) /
                                          1458. +
                                      (226. * (M_PI4)) / 405. + (2846. * z * (M_PI4)) / 1215.)
                                         .real();
    cache_ps[index_p(cu, 2, 2, 2)] = (27476329. / 58320. + (40370. * L_1) / 243. + (604. * L_12) / 27. + (6928. * L_2) / 81. +
                                      (1064. * L_1 * L_2) / 27. - (52. * L_22) / 3. + (416. * Cl2PI3 * sqrt3) / 81. -
                                      ((2. * complex::i()) / 27.) * log3 * log3 * sqrt3 - (26319. * sqrtz) / 250. - ((16. * -1.) / 27.) * sqrt3 * t_2 +
                                      (151070269. * z) / 81000. + (9466. * L_1 * z) / 27. + (464. * L_2 * z) / 9. + (464. * logz * z) / 3. -
                                      (4388. * zeta3) / 27. - 300. * z * zeta3 - (4. * sqrt3 * (20. + 30. * L_1 + 3. * log3 + 90. * z) * (M_PI)) / 243. - ((41279. - 1728. * log12sqrt52 + 11808. * log2 - (20. * complex::i()) * sqrt3 - 1872. * sqrt5 + 72342. * z - 7344. * log2 * z + 10368. * logz * z + 864. * sqrt5 * z + 11664. * L_1 * (1. + z) + 3456. * L_2 * (1. + z)) * (M_PI2)) / 486. + (398. * (M_PI4)) / 135. + (7991. * z * (M_PI4)) / 810.)
                                         .real();
 
    // Corrections from LO and NLO coefficients (incl. pengFlag)
    // resummation logz in NLO coefficients
    cache_p[index_p(cu, 1, 1, 2)] += ((logz * (33433. + 12132. * L_1 + 2376. * L_2 - 12. * M_PI2) - 22. * (151. + 9. * log2z + 12. * M_PI2)) * z) / 162.;
    cache_p[index_p(cu, 1, 2, 2)] += (-2. * (302. + 18. * log2z + logz * (1663. + 576. * L_1 - 216. * L_2 - 12. * M_PI2) + 24. * M_PI2) * z) / 27.;
    cache_p[index_p(cu, 2, 2, 2)] += (2. * (-151. - 9. * log2z + 2. * logz * (296. + 99. * L_1 + 54. * L_2 - 6. * M_PI2) - 12. * M_PI2) * z) / 9.;
    cache_ps[index_p(cu, 1, 1, 2)] += (4. * logz * (3473. + 12. * M_PI2) * z) / 81.;
    cache_ps[index_p(cu, 1, 2, 2)] += (-64. * logz * (-124. + 3. * M_PI2) * z) / 27.;
    cache_ps[index_p(cu, 2, 2, 2)] += (16. * logz * (-91. + 12. * M_PI2) * z) / 9.;
 
    // resummation affecting arguments of logs and Dilogs
    cache_p[index_p(cu, 1, 1, 2)] += 4762. / 81. * logz * z;
    cache_p[index_p(cu, 1, 2, 2)] += 1688. / 27. * logz * z;
    cache_p[index_p(cu, 2, 2, 2)] += 904. / 9. * logz * z;
    cache_ps[index_p(cu, 1, 1, 2)] += 16. / 81. * logz * z;
    cache_ps[index_p(cu, 1, 2, 2)] += -64. / 27. * logz * z;
    cache_ps[index_p(cu, 2, 2, 2)] += 64. / 9. * logz * z;
 
    // Transforming mb in NLO coefficients from Pole scheme to MSbar
    cache_p[index_p(cu, 1, 1, 2)] += (4. * (4. + 3. * L_b) * (-392. + 675. * z)) / 243.;
    cache_p[index_p(cu, 1, 2, 2)] += (-44. * (4. + 3. * L_b) * (-19. + 54. * z)) / 81.;
    cache_p[index_p(cu, 2, 2, 2)] += (-8. * (4. + 3. * L_b) * (-8. + 27. * z)) / 27.;
    cache_ps[index_p(cu, 1, 1, 2)] += (1264. * (4. + 3. * L_b)) / 243.;
    cache_ps[index_p(cu, 1, 2, 2)] += (416. * (4. + 3. * L_b)) / 81.;
    cache_ps[index_p(cu, 2, 2, 2)] += (-704. * (4. + 3. * L_b)) / 27.;
 
    for (quarks qq = cc; qq <= uu; qq = quarks(qq + 2))
    {
        // Gerlach thesis eq. (6.28)
        cache_p[index_p(qq, 1, 8, 2)] = 208. / 81. * L_1 - L_2 / 27. + (2615. / 54. - 10. / 9. * M_PI2) * z - 5. / 9. * M_PI2 + 25. / (54. * sqrt3) * M_PI - 115. / 486.;
        cache_p[index_p(qq, 2, 8, 2)] = -11. / 27. * L_1 + 2. / 9. * L_2 + (20. / 3. * M_PI2 - 833. / 9.) * z + 10. / 3. * M_PI2 - 25. / (9. * sqrt3) * M_PI - 3125. / 81.;
        cache_ps[index_p(qq, 1, 8, 2)] = 448. / 81. * L_1 + 32. / 27. * L_2 + (1192. / 27. - 16. / 9. * M_PI2) * z - 8. / 9. * M_PI2 + 20. / (27. * sqrt3) * M_PI + 3580. / 243.;
        cache_ps[index_p(qq, 2, 8, 2)] = -248. / 27. * L_1 - 64. / 9. * L_2 + (32. / 3. * M_PI2 - 1088. / 9.) * z + 16. / 3. * M_PI2 - 40. / (9. * sqrt3) * M_PI - 4568. / 81.;
 
        // equation (6.30)
        cache_p[index_p(qq, 3, 8, 2)] = -85. / 27. * L_1 - 448. / 9. * L_2 - 196. / 3. * z + 25. / 6. * M_PI2 - 107. / (18. * sqrt3) * M_PI - 17201. / 81.;
        cache_p[index_p(qq, 4, 8, 2)] = -3269. / 162. * L_1 - 427. / 27. * L_2 + (20. / 3. * M_PI2 - 404. / 3.) * z + 169. / 12. * M_PI2 - 514. / (27. * sqrt3) * M_PI - 43016. / 243.;
        cache_p[index_p(qq, 5, 8, 2)] = 5120. / 27. * L_1 - 7168. / 9. * L_2 - 760. / 3. * z + 770. / 9. * M_PI2 - 28. / (9. * sqrt3) * M_PI - 430238. / 81.;
        cache_p[index_p(qq, 6, 8, 2)] = -8962. / 81. * L_1 - 6976. / 27. * L_2 + (200. / 3. * M_PI2 - 4222. / 3.) * z + 3761. / 27. * M_PI2 - 3220. / (27. * sqrt3) * M_PI - 474656. / 243.;
        cache_ps[index_p(qq, 3, 8, 2)] = 440. / 27. * L_1 + 512. / 9. * L_2 + 608. / 3. * z - 596. / 27. * M_PI2 - 52. / (9. * sqrt3) * M_PI + 22504. / 81.;
        cache_ps[index_p(qq, 4, 8, 2)] = -4804. / 81. * L_1 - 160. / 27. * L_2 + (32. / 3. * M_PI2 - 128. / 3.) * z + 1090. / 81. * M_PI2 - 1120. / (27. * sqrt3) * M_PI - 46988. / 243.;
        cache_ps[index_p(qq, 5, 8, 2)] = 17408. / 27. * L_1 + 8192. / 9. * L_2 + 11456. / 3. * z - 8912. / 27. * M_PI2 - 1984. / (9. * sqrt3) * M_PI + 420304. / 81.;
        cache_ps[index_p(qq, 6, 8, 2)] = -28624. / 81. * L_1 + 2048. / 27. * L_2 + (320. / 3. * M_PI2 - 160. / 3.) * z + 5608. / 81. * M_PI2 - 8416. / (27. * sqrt3) * M_PI - 423440. / 243.;
 
        // Gerlach thesis eq. (6.31) and (6.29)
        z = 0.;
    }
    z = cache_z;
 
    // Gerlach thesis eq. (6.29)
    cache_p[index_p(uu, 1, 8, 2)] += -10. / 9. * z;
    cache_p[index_p(uu, 2, 8, 2)] += 20. / 3. * z;
    cache_ps[index_p(uu, 1, 8, 2)] += -16. / 9. * z;
    cache_ps[index_p(uu, 1, 8, 2)] += 32. / 3. * z;
 
    // Gerlach thesis eq. (6.31)
    cache_p[index_p(uu, 3, 8, 2)] += -196. / 3. * z;
    cache_p[index_p(uu, 4, 8, 2)] += (-404. / 3. + 20. / 3. * M_PI2) * z;
    cache_p[index_p(uu, 5, 8, 2)] += -760. / 3. * z;
    cache_p[index_p(uu, 6, 8, 2)] += (-4222. / 3. + 200. / 3. * M_PI2) * z;
    cache_ps[index_p(uu, 3, 8, 2)] += 608. / 3. * z;
    cache_ps[index_p(uu, 4, 8, 2)] += (-128. / 3. + 32. / 3. * M_PI2) * z;
    cache_ps[index_p(uu, 5, 8, 2)] += 11456. / 3. * z;
    cache_ps[index_p(uu, 6, 8, 2)] += (-160. / 3. + 320. / 3. * M_PI2) * z;
 
    // as in Gerlach thesis eq. (6.20)
    for (int i = 1; i <= 6; i++)
    {
        cache_p[index_p(cu, i, 8, 2)] =
            0.5 * (cache_p[index_p(cc, i, 8, 2)] + cache_p[index_p(uu, i, 8, 2)]);
        cache_ps[index_p(cu, i, 8, 2)] =
            0.5 * (cache_ps[index_p(cc, i, 8, 2)] + cache_ps[index_p(uu, i, 8, 2)]);
    }
 
    // Gerlach thesis eq. (6.32)
    cache_p[index_p(cc, 8, 8, 2)] = -13. / 18.;
    cache_p[index_p(cu, 8, 8, 2)] = -13. / 18.;
    cache_p[index_p(uu, 8, 8, 2)] = -13. / 18.;
    cache_ps[index_p(cc, 8, 8, 2)] = -68. / 9.;
    cache_ps[index_p(cu, 8, 8, 2)] = -68. / 9.;
    cache_ps[index_p(uu, 8, 8, 2)] = -68. / 9.;
 
    // compute flag_LOz terms
    for (int i = 0; i < 576; i++)
    {
        cache_p_LO[i] = cache_p[i];
        cache_ps_LO[i] = cache_ps[i];
    }
 
    cache_p_LO[index_p(cc, 1, 1, 0)] = 23. / 72. - 11. / 6. * z;
    cache_p_LO[index_p(cc, 1, 2, 0)] = 1. / 6. - 2. * z;
    cache_p_LO[index_p(cc, 2, 2, 0)] = 1. - 3. * z;
    cache_ps_LO[index_p(cc, 1, 1, 0)] = -5. / 9.;
    cache_ps_LO[index_p(cc, 1, 2, 0)] = -4. / 3.;
    cache_ps_LO[index_p(cc, 2, 2, 0)] = 1.;
 
    cache_p_LO[index_p(cc, 1, 5, 0)] = 64. / 3. - 96. * z;
    cache_p_LO[index_p(cc, 1, 6, 0)] = 4. * z - 20. / 9.;
    cache_p_LO[index_p(cc, 2, 5, 0)] = 16. - 72. * z;
    cache_p_LO[index_p(cc, 2, 6, 0)] = 40. / 3. - 24. * z;
 
    cache_p_LO[index_p(cc, 5, 5, 0)] = -1296. * n_v * z + 408. * (n_l + n_v) + 512.;
    cache_p_LO[index_p(cc, 6, 6, 0)] = -72. * n_v * z + 170. / 3. * (n_l + n_v) + 416. / 9.;
 
    for (int i = 3; i <= 6; i++)
    {
        for (int j = i; j <= 6; j++)
        {
            cache_p_LO[index_p(uu, i, j, 0)] = cache_p_LO[index_p(cc, i, j, 0)];
            cache_ps_LO[index_p(uu, i, j, 0)] = cache_ps_LO[index_p(cc, i, j, 0)];
        }
    }
 
    // Gerlach thesis eq. (6.10)
    for (int i = 1; i <= 6; i++)
    {
        for (int j = i; j <= 6; j++)
        {
            cache_p_LO[index_p(cu, i, j, 0)] =
                0.5 * (cache_p_LO[index_p(cc, i, j, 0)] + cache_p_LO[index_p(uu, i, j, 0)]);
            cache_ps_LO[index_p(cu, i, j, 0)] =
                0.5 * (cache_ps_LO[index_p(cc, i, j, 0)] + cache_ps_LO[index_p(uu, i, j, 0)]);
        }
    }
 
    cache_p_LO[index_p(cc, 1, 1, 1)] = (1196. + 342. * L_1 + 447. * L_2 - 15. * (M_PI2)) / 324. + (z * (-4113. - 1008. * L_1 - 792. * L_2 - 3168. * logz + 4. * (M_PI2))) / 216.;
    cache_p_LO[index_p(cc, 1, 2, 1)] = (z * (1179. + 468. * L_1 - 72. * L_2 - 288. * logz - 4. * (M_PI2))) / 18. + (-452. - (333. * L_1) / 2. + 57. * L_2 + 15. * (M_PI2)) / 27.;
    cache_p_LO[index_p(cc, 2, 2, 1)] = (37. - 18. * L_1 + 12. * L_2 - 30. * (M_PI2)) / 18. + (z * (135. + 72. * L_1 - 36. * L_2 - 144. * logz + 4. * (M_PI2))) / 6.;
    cache_ps_LO[index_p(cc, 1, 1, 1)] = z * (-385. / 9. - (4. * (M_PI2)) / 27.) - (2. * (-1. + 18. * L_1 + 60. * L_2 + 3. * (M_PI2))) / 81.;
    cache_ps_LO[index_p(cc, 1, 2, 1)] = (16. * z * (-42. + (M_PI2))) / 9. + (8. * (11. + (9. * L_1) / 2. - 12. * L_2 + 3. * (M_PI2))) / 27.;
    cache_ps_LO[index_p(cc, 2, 2, 1)] = z * (44. - (16. * (M_PI2)) / 3.) - (8. * (-31. - 9. * L_1 - 3. * L_2 + 3. * (M_PI2))) / 9.;
 
    cache_p_LO[index_p(cc, 1, 1, 1)] += -5. / 486. - (5. * L_1) / 324. - (5. * z) / 54.;
    cache_p_LO[index_p(cc, 1, 2, 1)] += 10. / 81. + (5. * L_1) / 27. + (10. * z) / 9.;
    cache_p_LO[index_p(cc, 2, 2, 1)] += -10. / 27. - (5. * L_1) / 9. - (10. * z) / 3.;
    cache_ps_LO[index_p(cc, 1, 1, 1)] += -4. / 243. - (2. * L_1) / 81. - (4. * z) / 27.;
    cache_ps_LO[index_p(cc, 1, 2, 1)] += 16. / 81. + (8. * L_1) / 27. + (16. * z) / 9.;
    cache_ps_LO[index_p(cc, 2, 2, 1)] += -16. / 27. - (8. * L_1) / 9. - (16. * z) / 3.;
 
    cache_p_LO[index_p(cc, 1, 1, 1)] += -11. / 6. * z_replace;
    cache_p_LO[index_p(cc, 1, 2, 1)] += -2. * z_replace;
    cache_p_LO[index_p(cc, 2, 2, 1)] += -3. * z_replace;
 
    // Gerlach thesis eq. (6.13) and (6.16)
    for (int i = 1; i <= 2; i++)
    {
        for (int j = i; j <= 2; j++)
        {
            cache_p_LO[index_p(cu, i, j, 1)] =
                0.5 * (cache_p_LO[index_p(cc, i, j, 1)] + cache_p_LO[index_p(uu, i, j, 1)]);
            cache_ps_LO[index_p(cu, i, j, 1)] =
                0.5 * (cache_ps_LO[index_p(cc, i, j, 1)] + cache_ps_LO[index_p(uu, i, j, 1)]);
        }
    }
 
    // Gerlach thesis eq. (6.19)
    cache_p_LO[index_p(cc, 1, 8, 1)] = 5. / 18.;
    cache_p_LO[index_p(cc, 2, 8, 1)] = -5. / 3.;
    cache_ps_LO[index_p(cc, 1, 8, 1)] = 4. / 9.;
    cache_ps_LO[index_p(cc, 2, 8, 1)] = -8. / 3.;
 
    // Gerlach thesis eq. (6.21)
    cache_p_LO[index_p(cc, 3, 8, 1)] = -32. / 3.;
    cache_p_LO[index_p(cc, 4, 8, 1)] = -169. / 18.;
    cache_p_LO[index_p(cc, 5, 8, 1)] = -512. / 3.;
    cache_p_LO[index_p(cc, 6, 8, 1)] = -992. / 9.;
    cache_ps_LO[index_p(cc, 3, 8, 1)] = 64. / 3.;
    cache_ps_LO[index_p(cc, 4, 8, 1)] = -20. / 9.;
    cache_ps_LO[index_p(cc, 5, 8, 1)] = 1024. / 3.;
    cache_ps_LO[index_p(cc, 6, 8, 1)] = 256. / 9.;
 
    // Gerlach thesis eq. (6.20)
    for (int i = 1; i <= 6; i++)
    {
        cache_p_LO[index_p(cu, i, 8, 1)] = cache_p_LO[index_p(cc, i, 8, 1)];
        cache_p_LO[index_p(uu, i, 8, 1)] = cache_p_LO[index_p(cc, i, 8, 1)];
        cache_ps_LO[index_p(cu, i, 8, 1)] = cache_ps_LO[index_p(cc, i, 8, 1)];
        cache_ps_LO[index_p(uu, i, 8, 1)] = cache_ps_LO[index_p(cc, i, 8, 1)];
    }
 
    logz = cache_logz;
    //lastInput_compute_pp_s[0] = z;
    //lastInput_compute_pp_s[1] = mu_1;
    //lastInput_compute_pp_s[2] = mu_2;
    //lastInput_compute_pp_s[3] = mu_b;
    return;
}

◆ computeCKMandMasses()

void AmpDB2::computeCKMandMasses	(	orders	order = `NNLO`,
		mass_schemes	mass_scheme = `MSbar`
	)

private

A method to compute CKM elements, quark masses and alpha_s.

Parameters

[in]	order	NNLO, use NLO only for tradBasis
[in]	mass_scheme	the scheme for the bottom quark mass

Definition at line 317 of file AmpDB2.cpp.

{
    if (order != NLO and order != NNLO)
        throw(std::runtime_error("computeCKMandMasses() order not present"));
 
    VtbVtd = mySM.getCKM().computelamt_d();
    VtbVts = mySM.getCKM().computelamt_s();
    VtbVtd2 = VtbVtd * VtbVtd;
    VtbVts2 = VtbVts * VtbVts;
    VcbVcd = mySM.getCKM().computelamc_d();
    VcbVcs = mySM.getCKM().computelamc_s();
    VcbVcd2 = VcbVcd * VcbVcd;
    VcbVcs2 = VcbVcs * VcbVcs;
    lambda_c_d = VcbVcd.conjugate();
    lambda_u_d = mySM.getCKM().computelamu_d().conjugate();
    lambda_c_s = VcbVcs.conjugate();
    lambda_u_s = mySM.getCKM().computelamu_s().conjugate();
 
    // pole mass of bottom quark
    Mb_pole = mySM.Mbar2Mp(mySM.getQuarks(QCD::BOTTOM).getMass(), QCD::BOTTOM);
 
    // DB=1 matching scales (arxiv: 2205.07907 Results. or Gerlach thesis eq. 7.7) varied by "getMub()" or fixed to 4.2
    mu_1 = mySM.getMub();
    // mass scale mu_b=mu_c
    mu_b = mu_1;
    if (flag_fixmub)
        mu_b = 4.2;
    // mu_1_overm = Mb_pole; // Gerlach thesis
    mu_1_overm = mu_1;
 
    // DB=2 matching scale mu_2
    mu_2 = mySM.getBBs().getMu();
 
    // MSbar bottom quark mass Mb(Mb)
    Mb_Mb = mySM.getQuarks(QCD::BOTTOM).getMass();
 
    // MSbar bottom quark mass Mb(mu_b)
    double Mb_mub = mySM.Mrun(mu_b,
                              mySM.getQuarks(QCD::BOTTOM).getMass_scale(),
                              mySM.getQuarks(QCD::BOTTOM).getMass(), QCD::BOTTOM, FULLNNLO);
 
    // MSbar charm quark mass Mc(Mc)
    double Mc_Mc = mySM.getQuarks(QCD::CHARM).getMass();
 
    // MSbar charm quark mass Mc(mu_b)
    double Mc_mub = mySM.Mrun(mu_b,
                              mySM.getQuarks(QCD::CHARM).getMass_scale(),
                              mySM.getQuarks(QCD::CHARM).getMass(), QCD::CHARM, FULLNNLO);
 
    // strong coupling constant divided by 4*Pi
    as_4pi_mu1 = mySM.Als(mu_1, FULLNNNLO, true, true) / (4. * M_PI);
    as_4pi_mu2 = mySM.Als(mu_2, FULLNNNLO, true, true) / (4. * M_PI);
    as_4pi = mySM.Als(Mb_Mb, FULLNNNLO, true, true) / (4. * M_PI);
 
    // resummed z always in MSbar
    z = Mc_mub * Mc_mub / (Mb_mub * Mb_mub);
    sqrtz = sqrt(z);
 
    // PS mass of bottom quark: NNLO evaluation from hep-ph/9804241v2 eq. (25)
    PoletoMS_as2 = 307. / 32. + M_PI2 / 3. + M_PI2 / 9. * log2 - 1. / 6. * zeta3 - (71. / 144. + M_PI2 / 18.) * nl + 4. / 3. * M_PI2 / 8. * sqrtz - z;
    Mb_PS = Mb_Mb * (1. + 16. / 3. * as_4pi * (1. - mu_f / Mb_Mb) + 16. * as_4pi * as_4pi * (PoletoMS_as2 - mu_f / (3. * Mb_Mb) * (a1 - b0 * (2. * log(mu_f / Mb_Mb) - 2.))));
 
    // adapt "Mb2_prefactor" to the used mass scheme; Mb, Mc, always in MSbar
    // explained in Gerlach thesis chapter 7.0 and arxiv:2205.07907 Results.
    if (order == NNLO)
    {
        Mb = Mb_mub; // if Mb changes independently of mub it has to be added to currentInput_compute_pp_s
        Mc = Mc_mub;
        this->flag_resumz = true;
        switch (mass_scheme)
        {
        case pole:
            Mb2_prefactor = Mb_pole * Mb_pole;
            break;
        case MSbar_partialNNLO:
        case MSbar_partialN3LO:
        case MSbar:
            Mb2_prefactor = Mb_mub * Mb_mub;
            break;
        case PS_partialNNLO:
        case PS_partialN3LO:
        case PS:
            Mb2_prefactor = Mb_PS * Mb_PS;
            break;
        case only1overmb:
            Mb2_prefactor = 0.;
            break;
        default:
            throw(std::runtime_error("mass_scheme not implemented"));
        }
    }
    Mb2_prefactor_1overm = Mb2_prefactor;
    if (mass_scheme == only1overmb) 
        Mb2_prefactor_1overm = Mb_PS * Mb_PS;
 
    // adapt to MSbar mass scheme for the traditional basis
    if (order == NLO)
    {
        Mb2_prefactor = Mb_mub * Mb_mub;
        Mb = Mb_pole;
        x_1 = mu_1 / Mb;
        x_2 = mu_2 / Mb;
        logx_1 = log(x_1);
        logx_2 = log(x_2);
        double Mc = Mc_mub;
        if (!flag_resumz)
            Mc = Mc_Mc;
        z = Mc * Mc / (Mb_mub * Mb_mub);
    }
 
    Gf2 = mySM.getGF() * mySM.getGF();
    MB = mySM.getMesons(QCD::B_D).getMass();
    MB_s = mySM.getMesons(QCD::B_S).getMass();
 
    z2 = z * z;
    logz = log(z);
    log2z = logz * logz;
    oneminusz2 = (1. - z) * (1. - z);
    sqrt1minus4z = sqrt(1. - 4. * z);
 
    // calculate "z" values for 1/m_b contributions
    // z_1overm = Mc_Mc * Mc_Mc / (Mb_Mb * Mb_Mb); //hep-ph/0612167 eq. 27
    // MSbar charm quark mass Mb(Mb)
    // double Mc_Mb = mySM.Mrun(mySM.getQuarks(QCD::BOTTOM).getMass_scale(),
    //     mySM.getQuarks(QCD::CHARM).getMass_scale(),
    //     mySM.getQuarks(QCD::CHARM).getMass(), FULLNNLO);
    // z_1overm = Mc_Mb * Mc_Mb / (Mb_Mb * Mb_Mb); //arxiv:1910.00970 eq. 11
    // if z for the 1/m_b contributions shall be varied with "mu_b"
    z_1overm = Mc_mub * Mc_mub / (Mb_mub * Mb_mub);
 
    z_1overm2 = z_1overm * z_1overm;
    oneminusz_1overm2 = (1. - z_1overm) * (1. - z_1overm);
    sqrt1minus4z_1overm = sqrt(1. - 4. * z_1overm);
 
    z3 = z2 * z;
    z4 = z3 * z;
 
    log1minusz = log(1. - z);
    log1minus4z = log(1. - 4. * z);
 
    sigma = (1. - sqrt1minus4z) / (1. + sqrt1minus4z);
    logsigma = log(sigma);
    log2sigma = logsigma * logsigma;
 
    Dilogz = gslpp_special_functions::dilog(z);
    Dilogsigma = gslpp_special_functions::dilog(sigma);
    Dilogsigma2 = gslpp_special_functions::dilog(sigma * sigma);
    return;
}

◆ computeD()

void AmpDB2::computeD ( orders order )

private

Definition at line 640 of file AmpDB2.cpp.

{
    if (order == FULLNLO)
    {
        // qq = uu and cc
        for (quarks qq = cc; qq <= uu; qq = quarks(qq + 2))
        {
            for (int k = 1; k <= 2; k++)
            {
                gslpp::complex result = 0.;
                for (int i = 1; i <= 2; i++)
                {
                    for (int j = 1; j <= 2; j++)
                    {
                        result += C_Buras_LO[i - 1] * C_Buras_LO[j - 1] * F(qq, k, i, j) + (C_Buras_NLO[i - 1] * C_Buras_LO[j - 1] + C_Buras_LO[i - 1] * C_Buras_NLO[j - 1]) * F0(qq, k, i, j);
                    }
                }
                result += +as_4pi_mu1 * C_Buras_LO[2 - 1] * C_Buras_LO[2 - 1] * P(qq, k, 2, 2) + 2. * as_4pi_mu1 * C_Buras_LO[2 - 1] * C_Buras_LO[8 - 1] * P(qq, k, 2, 8);
                for (int i = 1; i <= 2; i++)
                {
                    for (int r = 3; r <= 6; r++)
                    {
                        result += 2. * C_Buras_LO[i - 1] * C_Buras_LO[r - 1] * P(qq, k, i, r);
                    }
                }
                cacheD[indexD(qq, k)] = result;
            }
        }
        // qq = cu
        for (int k = 1; k <= 2; k++)
        {
            gslpp::complex result = 0.;
            for (int i = 1; i <= 2; i++)
            {
                for (int j = 1; j <= 2; j++)
                {
                    result += C_Buras_LO[i - 1] * C_Buras_LO[j - 1] * F(cu, k, i, j) + (C_Buras_NLO[i - 1] * C_Buras_LO[j - 1] + C_Buras_LO[i - 1] * C_Buras_NLO[j - 1]) * F0(cu, k, i, j);
                }
            }
            result += +as_4pi_mu1 * C_Buras_LO[2 - 1] * C_Buras_LO[2 - 1] * P(cu, k, 2, 2) + as_4pi_mu1 * C_Buras_LO[2 - 1] * C_Buras_LO[8 - 1] * (P(cc, k, 2, 8) + P(uu, k, 2, 8));
            for (int i = 1; i <= 2; i++)
            {
                for (int r = 3; r <= 6; r++)
                {
                    result += C_Buras_LO[i - 1] * C_Buras_LO[r - 1] * (P(cc, k, i, r) + P(uu, k, i, r));
                }
            }
            cacheD[indexD(cu, k)] = result;
        }
    }
    else if (order == LO)
    {
        // qq = uu and cc
        for (quarks qq = cc; qq <= uu; qq = quarks(qq + 2))
        {
            for (int k = 1; k <= 2; k++)
            {
                gslpp::complex result = 0.;
                for (int i = 1; i <= 2; i++)
                {
                    for (int j = 1; j <= 2; j++)
                    {
                        result += C_Buras_LO[i - 1] * C_Buras_LO[j - 1] * F0(qq, k, i, j);
                    }
                }
                for (int i = 1; i <= 2; i++)
                {
                    for (int r = 3; r <= 6; r++)
                    {
                        result += 2. * C_Buras_LO[i - 1] * C_Buras_LO[r - 1] * P(qq, k, i, r);
                    }
                }
                cacheD[indexD(qq, k)] = result;
            }
        }
        // qq = cu
        for (int k = 1; k <= 2; k++)
        {
            gslpp::complex result = 0.;
            for (int i = 1; i <= 2; i++)
            {
                for (int j = 1; j <= 2; j++)
                {
                    result += C_Buras_LO[i - 1] * C_Buras_LO[j - 1] * F0(cu, k, i, j);
                }
            }
            for (int i = 1; i <= 2; i++)
            {
                for (int r = 3; r <= 6; r++)
                {
                    result += C_Buras_LO[i - 1] * C_Buras_LO[r - 1] * (P(cc, k, i, r) + P(uu, k, i, r));
                }
            }
            cacheD[indexD(cu, k)] = result;
        }
    }
    else
    {
        throw std::runtime_error("AmpDB2::computeD(orders order): order not implemented");
    }
    return;
}

◆ computeF0()

void AmpDB2::computeF0 ( )

private

Definition at line 490 of file AmpDB2.cpp.

{
    cacheF0[indexF(cu, 1, 1, 1)] = 1.5 * (2. - 3. * z + z3);
    cacheF0[indexF(cu, 1, 1, 2)] = 0.5 * (2. - 3. * z + z3);
    cacheF0[indexF(cu, 1, 2, 2)] = 0.5 * oneminusz2 * (1. - z);
    cacheF0[indexF(cu, 2, 1, 1)] = 3. * oneminusz2 * (1. + 2. * z);
    cacheF0[indexF(cu, 2, 1, 2)] = oneminusz2 * (1. + 2. * z);
    cacheF0[indexF(cu, 2, 2, 2)] = -oneminusz2 * (1. + 2. * z);
    cacheF0[indexF(cc, 1, 1, 1)] = 3. * sqrt1minus4z * (1. - z);
    cacheF0[indexF(cc, 1, 1, 2)] = sqrt1minus4z * (1. - z);
    cacheF0[indexF(cc, 1, 2, 2)] = 0.5 * sqrt1minus4z * sqrt1minus4z * sqrt1minus4z;
    cacheF0[indexF(cc, 2, 1, 1)] = 3. * sqrt1minus4z * (1. + 2. * z);
    cacheF0[indexF(cc, 2, 1, 2)] = sqrt1minus4z * (1. + 2. * z);
    cacheF0[indexF(cc, 2, 2, 2)] = -sqrt1minus4z * (1. + 2. * z);
    cacheF0[indexF(uu, 1, 1, 1)] = 3.;
    cacheF0[indexF(uu, 1, 1, 2)] = 1.;
    cacheF0[indexF(uu, 1, 2, 2)] = 0.5;
    cacheF0[indexF(uu, 2, 1, 1)] = 3.;
    cacheF0[indexF(uu, 2, 1, 2)] = 1.;
    cacheF0[indexF(uu, 2, 2, 2)] = -1.;
    for (int k = 1; k < 3; k++)
    {
        for (quarks qq = cc; qq <= uu; qq = quarks(qq + 1))
        {
            cacheF0[indexF(qq, k, 2, 1)] = cacheF0[indexF(qq, k, 1, 2)];
        }
    }
    return;
}

◆ computeF1()

void AmpDB2::computeF1 ( )

private

Definition at line 523 of file AmpDB2.cpp.

{
    double log_muM = 2. * log(mu_b / Mb_pole);
    cacheF1[indexF(cu, 1, 1, 1)] = 109. / 6. - 37. * z + 1.5 * z2 + 52. / 3. * z3 + 2. * oneminusz2 * (5. + z) * logx_2 - 4. * oneminusz2 * (5. + 7. * z) * log1minusz -
                                   2. * z * (10. + 14. * z - 15. * z2) * logz + 8. * (2. - 3. * z + z3) * log1minusz * logz + 16. * (2. - 3. * z + z3) * Dilogz + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(cu, 1, 1, 1) - 8. * z * logz * 1.5 * (-3. + 3. * z2));
    cacheF1[indexF(cu, 2, 1, 1)] = -4. / 3. * (10. - 33. * z + 54. * z2 - 31. * z3) - 8. * oneminusz2 * (4. + 14. * z - 3. * z2) * log1minusz +
                                   8. * z * (2. - 23. * z + 21. * z2 - 3. * z3) * logz -
                                   16. * oneminusz2 * (1. + 2. * z) * (2. * logx_2 - log1minusz * logz - 2. * Dilogz) + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(cu, 2, 1, 1) - 8. * z * logz * 3. * 6. * (z - 1.) * z);
    cacheF1[indexF(cu, 1, 1, 2)] = 1. * ((-502. + 912. * z - 387. * z2 - 23. * z3) / 36. - oneminusz2 * (17. + 4. * z) * logx_1 + 2. / 3. * oneminusz2 * (5. + z) * logx_2 -
                                         oneminusz2 / (12. * z) * (2. + 33. * z + 94. * z2) * log1minusz - z / 12. * (80. + 69. * z - 126. * z2) * logz +
                                         8. / 3. * (2. - 3. * z + z3) * (log1minusz * logz + 2. * Dilogz)) +
                                   flag_resumz * ((32. / 3. + 8. * log_muM) * F0(cu, 1, 1, 2) - 8. * z * logz * 0.5 * (-3. + 3. * z2));
    cacheF1[indexF(cu, 2, 1, 2)] = 1. * ((-130. + 93. * z + 144. * z2 - 107. * z3) / 9. - 2. / 3. * oneminusz2 / z * (1. + 15. * z + 47. * z2 - 12. * z3) * log1minusz +
                                         2. / 3. * z * (8. - 93. * z + 87. * z2 - 12. * z3) * logz -
                                         8. / 3. * oneminusz2 * (1. + 2. * z) * (3. * logx_1 + 4. * logx_2 - 2. * log1minusz * logz - 4. * Dilogz)) +
                                   flag_resumz * ((32. / 3. + 8. * log_muM) * F0(cu, 2, 1, 2) - 8. * z * logz * 6. * (z - 1.) * z);
    cacheF1[indexF(cu, 1, 2, 2)] = -M_PI2 / 3. * (1. - 5. * z + 4. * z2) + (-136. - 159. * z + 738. * z2 - 443. * z3) / 18. - 2 * oneminusz2 * (5. + 4. * z) * logx_1 +
                                   2. / 3. * oneminusz2 * (4. - z) * logx_2 + oneminusz2 / (6. * z) * (7. + 32. * z2 + 3. * z3) * log1minusz -
                                   z / 6. * (62. + 39. * z - 30. * z2 + 3. * z3) * logz + (5. - 3. * z - 18. * z2 + 16. * z3) / 3. * (log1minusz * logz + 2. * Dilogz) + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(cu, 1, 2, 2) - 8. * z * logz * -1.5 * oneminusz2);
    cacheF1[indexF(cu, 2, 2, 2)] = 8. / 3. * M_PI2 * (1. + z - 2. * z2) - 28. / 9. * (5. + 3. * z - 27. * z2 + 19. * z3) - 16. * oneminusz2 * (1. + 2. * z) * logx_1 +
                                   32. / 3. * oneminusz2 * (1. + 2. * z) * logx_2 - 4. / 3. * oneminusz2 / z * (1. - 12. * z - 16. * z2 - 3. * z3) * log1minusz +
                                   4. / 3. * z * (2. - 3. * z + 18. * z2 - 3. * z3) * logz + 8. / 3. * (1. - 3. * z - 6. * z2 + 8. * z3) * (log1minusz * logz + 2. * Dilogz) + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(cu, 2, 2, 2) - 8. * z * logz * -6. * (z - 1.) * z);
 
    cacheF1[indexF(cc, 1, 1, 1)] = sqrt1minus4z * (109. - 226. * z + 168. * z2) / 6. - (52. - 104. * z - 16. * z2 + 56. * z3) * logsigma +
                                   2. * (5. - 8. * z) * sqrt1minus4z * logx_2 - 12. * sqrt1minus4z * (3. - 2. * z) * log1minus4z + 4. * (13. - 10. * z) * sqrt1minus4z * logz +
                                   16. * (1. - 3. * z + 2. * z2) * (3. * log2sigma + 2. * logsigma * log1minus4z - 3. * logsigma * logz + 4. * Dilogsigma + 2. * Dilogsigma2) + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(cc, 1, 1, 1) - 8. * z * logz * 3. * (6. * z - 3.) / sqrt1minus4z);
    cacheF1[indexF(cc, 2, 1, 1)] = -8. / 3. * sqrt1minus4z * (5. - 23. * z - 42. * z2) - 16. * (4. - 2. * z - 7. * z2 + 14. * z3) * logsigma - 32. * sqrt1minus4z * (1. + 2. * z) * logx_2 -
                                   48. * sqrt1minus4z * (1. + 2. * z) * log1minus4z + 64. * sqrt1minus4z * (1. + 2. * z) * logz +
                                   16. * (1. - 4. * z2) * (3. * log2sigma + 2. * logsigma * log1minus4z - 3. * logsigma * logz + 4. * Dilogsigma + 2. * Dilogsigma2) + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(cc, 2, 1, 1) - 8. * z * logz * 3. * -12. * z / sqrt1minus4z);
    cacheF1[indexF(cc, 1, 1, 2)] = -sqrt1minus4z * (127. - 199. * z - 75. * z2) / 9. + (2. - 259. * z + 662. * z2 - 76. * z3 - 200. * z4) * logsigma / (12. * z) -
                                   (17. - 26. * z) * sqrt1minus4z * logx_1 + 2. / 3. * (5. - 8. * z) * sqrt1minus4z * logx_2 - 4. * sqrt1minus4z * (3. - 2. * z) * log1minus4z -
                                   sqrt1minus4z * (2. - 255. * z + 316. * z2) * logz / (12. * z) +
                                   16. / 3. * (1. - 3. * z + 2. * z2) * (3. * log2sigma + 2. * logsigma * log1minus4z - 3. * logsigma * logz + 4. * Dilogsigma + 2. * Dilogsigma2) + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(cc, 1, 1, 2) - 8. * z * logz * (6. * z - 3.) / sqrt1minus4z);
    cacheF1[indexF(cc, 2, 1, 2)] = -2. * sqrt1minus4z * (68. + 49. * z - 150. * z2) / 9. + 2. / 3. * (1. - 35. * z + 4. * z2 + 76. * z3 - 100. * z4) * logsigma / z +
                                   (16. - 64. * z2) * log2sigma - 8. * sqrt1minus4z * (1. + 2. * z) * logx_1 - 32. / 3. * sqrt1minus4z * (1. + 2. * z) * logx_2 -
                                   16. * sqrt1minus4z * (1. + 2. * z) * log1minus4z - 2. / 3. * sqrt1minus4z * (1. - 33. * z - 76. * z2) * logz / z +
                                   16. / 3. * (1. - 4. * z2) * (2. * logsigma * log1minus4z - 3. * logsigma * logz + 4. * Dilogsigma + 2. * Dilogsigma2) + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(cc, 2, 1, 2) - 8. * z * logz * -12. * z / sqrt1minus4z);
    cacheF1[indexF(cc, 1, 2, 2)] = -M_PI2 / 3. * (1. - 10. * z) - sqrt1minus4z * (115. + 632. * z + 96. * z2) / 18. - (7. + 13. * z - 194. * z2 + 304. * z3 - 64. * z4) * logsigma / (6. * z) -
                                   2. * sqrt1minus4z * (5. - 2. * z) * logx_1 + 4. / 3. * (2. - 5. * z) * sqrt1minus4z * logx_2 - 4. * (1. - 6. * z) * sqrt1minus4z * log1minus4z +
                                   (13. - 54. * z + 8. * z2) * logsigma * log1minus4z / 3. + sqrt1minus4z * (7. + 27. * z - 250. * z2) * logz / (6. * z) +
                                   (7. - 32. * z + 4. * z2) * (log2sigma - logsigma * logz) + 4. / 3. * (5. - 12. * z + 4. * z2) * Dilogsigma + 4. / 3. * (4. - 21. * z + 2. * z2) * Dilogsigma2 + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(cc, 1, 2, 2) - 8. * z * logz * 0.5 * -6. * sqrt1minus4z);
    cacheF1[indexF(cc, 2, 2, 2)] = 8. / 3. * M_PI2 * (1. + 2. * z) - 8. / 9. * sqrt1minus4z * (19. + 53. * z + 24. * z2) + 4. / 3. * (1. + 7. * z + 10. * z2 - 68. * z3 + 32. * z4) * logsigma / z -
                                   8. * (1. + 2. * z) * (1. + 2. * z) * log2sigma - 16. * sqrt1minus4z * (1. + 2. * z) * logx_1 + 32. / 3. * sqrt1minus4z * (1. + 2. * z) * logx_2 +
                                   16. * sqrt1minus4z * (1. + 2. * z) * log1minus4z - 8. / 3. * (1. + 6. * z + 8. * z2) * logsigma * log1minus4z -
                                   4. / 3. * sqrt1minus4z * (1. + 9. * z + 26. * z2) * logz / z + 8. * (1. + 2. * z) * (1. + 2. * z) * logsigma * logz +
                                   32. / 3. * (1. - 4. * z2) * Dilogsigma - 32. / 3. * (1. + 3. * z + 2. * z2) * Dilogsigma2 + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(cc, 2, 2, 2) - 8. * z * logz * 12. * z / sqrt1minus4z);
 
    cacheF1[indexF(uu, 1, 1, 1)] = 109. / 6. + 10. * logx_2 + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(uu, 1, 1, 1));
    cacheF1[indexF(uu, 2, 1, 1)] = -40. / 3. - 32. * logx_2 + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(uu, 2, 1, 1));
    cacheF1[indexF(uu, 1, 1, 2)] = -127. / 9. + 4. / 12. - 17. * logx_1 + 10. / 3. * logx_2 + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(uu, 1, 1, 2));
    cacheF1[indexF(uu, 2, 1, 2)] = -(136. + 12.) / 9. - 8. * logx_1 - 32. / 3. * logx_2 + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(uu, 2, 1, 2));
    cacheF1[indexF(uu, 1, 2, 2)] = -M_PI2 / 3. - (115. + 42.) / 18. - 10. * logx_1 + 8. / 3. * logx_2 + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(uu, 1, 2, 2));
    cacheF1[indexF(uu, 2, 2, 2)] = 8. * M_PI2 / 3. - 8. / 9. * (19. - 3.) - 16. * logx_1 + 32. / 3. * logx_2 + flag_resumz * ((32. / 3. + 8. * log_muM) * F0(uu, 2, 2, 2));
 
    for (int k = 1; k < 3; k++)
    {
        for (quarks qq = cc; qq <= uu; qq = quarks(qq + 1))
        {
            cacheF1[indexF(qq, k, 2, 1)] = cacheF1[indexF(qq, k, 1, 2)];
        }
    }
    return;
}

◆ computeP()

void AmpDB2::computeP ( )

private

Definition at line 588 of file AmpDB2.cpp.

{
    cacheP[indexP(cu, 1, 2, 2)] = -1. / 27. + 2. / 9. * z - logx_1 / 9. - sqrt1minus4z * (1. + 2. * z) * (2. + 3. * logsigma + 6. * logx_1) / 54. + logz / 18.;
    cacheP[indexP(cu, 2, 2, 2)] = 8. / 27. + 16. / 9. * z + 8. / 9. * logx_1 + 4. / 27. * sqrt1minus4z * (1. + 2. * z) * (2. + 3. * logsigma + 6. * logx_1) - 4. * logz / 9.;
    cacheP[indexP(cc, 1, 2, 2)] = -2. / 27. * sqrt1minus4z * (1. + 8. * z + 12. * z2) - logsigma / 9. + 4. / 3. * z2 * logsigma + 16. / 9. * z3 * logsigma -
                                  sqrt1minus4z * (1. + 2. * z) * (2. * logx_1 - logz) / 9.;
    cacheP[indexP(cc, 2, 2, 2)] = 16. / 27. * sqrt1minus4z * (1. + 8. * z + 12. * z2) + 8. / 9. * logsigma - 32. / 3. * z2 * logsigma - 128. / 9. * z3 * logsigma +
                                  8. / 9. * sqrt1minus4z * (1. + 2. * z) * (2. * logx_1 - logz);
    cacheP[indexP(uu, 1, 2, 2)] = -2. / 27. - 2. / 9. * logx_1;
    cacheP[indexP(uu, 2, 2, 2)] = 16. / 27. + 16. / 9. * logx_1;
 
    cacheP[indexP(cc, 1, 1, 3)] = 3. * sqrt1minus4z * (1. - z);
    cacheP[indexP(cc, 2, 1, 3)] = 3. * sqrt1minus4z * (1. + 2. * z);
    cacheP[indexP(cc, 1, 2, 3)] = sqrt1minus4z * (1. - z);
    cacheP[indexP(cc, 2, 2, 3)] = sqrt1minus4z * (1. + 2. * z);
    cacheP[indexP(cc, 1, 1, 4)] = sqrt1minus4z * (1. - z);
    cacheP[indexP(cc, 2, 1, 4)] = sqrt1minus4z * (1. + 2. * z);
    cacheP[indexP(cc, 1, 2, 4)] = 0.5 * sqrt1minus4z * sqrt1minus4z * sqrt1minus4z;
    cacheP[indexP(cc, 2, 2, 4)] = -sqrt1minus4z * (1. + 2. * z);
    cacheP[indexP(cc, 1, 1, 5)] = 9. * z * sqrt1minus4z;
    cacheP[indexP(cc, 2, 1, 5)] = 0.;
    cacheP[indexP(cc, 1, 2, 5)] = 3. * z * sqrt1minus4z;
    cacheP[indexP(cc, 2, 2, 5)] = 0.;
    cacheP[indexP(cc, 1, 1, 6)] = 3. * z * sqrt1minus4z;
    cacheP[indexP(cc, 2, 1, 6)] = 0.;
    cacheP[indexP(cc, 1, 2, 6)] = 3. * z * sqrt1minus4z;
    cacheP[indexP(cc, 2, 2, 6)] = 0.;
    cacheP[indexP(cc, 1, 2, 8)] = -1. / 6. * sqrt1minus4z * (1. + 2. * z);
    cacheP[indexP(cc, 2, 2, 8)] = 4. / 3. * sqrt1minus4z * (1. + 2. * z);
 
    cacheP[indexP(uu, 1, 1, 3)] = 3.;
    cacheP[indexP(uu, 2, 1, 3)] = 3.;
    cacheP[indexP(uu, 1, 2, 3)] = 1.;
    cacheP[indexP(uu, 2, 2, 3)] = 1.;
    cacheP[indexP(uu, 1, 1, 4)] = 1.;
    cacheP[indexP(uu, 2, 1, 4)] = 1.;
    cacheP[indexP(uu, 1, 2, 4)] = 0.5;
    cacheP[indexP(uu, 2, 2, 4)] = -1.;
    cacheP[indexP(uu, 1, 1, 5)] = 0.;
    cacheP[indexP(uu, 2, 1, 5)] = 0.;
    cacheP[indexP(uu, 1, 2, 5)] = 0.;
    cacheP[indexP(uu, 2, 2, 5)] = 0.;
    cacheP[indexP(uu, 1, 1, 6)] = 0.;
    cacheP[indexP(uu, 2, 1, 6)] = 0.;
    cacheP[indexP(uu, 1, 2, 6)] = 0.;
    cacheP[indexP(uu, 2, 2, 6)] = 0.;
    cacheP[indexP(uu, 1, 2, 8)] = -1. / 6.;
    cacheP[indexP(uu, 2, 2, 8)] = 4. / 3.;
    return;
}

◆ computeWilsonCoeffsBuras()

void AmpDB2::computeWilsonCoeffsBuras ( )

private

Definition at line 467 of file AmpDB2.cpp.

{
    // NNLO DB=1 Wilson coefficients
    gslpp::vector<gslpp::complex> **WilsonCoeffsBuras = mySM.getFlavour().ComputeCoeffsgamma_Buras(mu_1);
    for (int i = 0; i < 8; i++)
    {
        if (i == 6)
            i = 7;
        C_Buras_LO[i] = (*(WilsonCoeffsBuras[LO]))(i);
        C_Buras_NLO[i] = (*(WilsonCoeffsBuras[NLO]))(i);
        C_Buras_NNLO[i] = (*(WilsonCoeffsBuras[NNLO]))(i);
    }
 
    // LO DB=1 Wilson coefficients for 1/mb corrections
    gslpp::vector<gslpp::complex> **WilsonCoeffs_1overm = mySM.getFlavour().ComputeCoeffsgamma_Buras(mu_1_overm);
    C1_LO_1overm = (*(WilsonCoeffs_1overm[LO]))(0);
    C2_LO_1overm = (*(WilsonCoeffs_1overm[LO]))(1);
    K_1 = 3. * C1_LO_1overm * C1_LO_1overm + 2. * C1_LO_1overm * C2_LO_1overm;
    K_2 = C2_LO_1overm * C2_LO_1overm;
    return;
}

◆ computeWilsonCoeffsMisiak()

void AmpDB2::computeWilsonCoeffsMisiak ( )

private

Definition at line 1173 of file AmpDB2.cpp.

{
    // NLO DB=1 Wilson coefficients C_i, i=1-6,8
    gslpp::vector<gslpp::complex> **WilsonCoeffsMisiak = mySM.getFlavour().ComputeCoeffsgamma(mu_1);
    for (int i = 0; i < 8; i++)
    {
        if (i == 6)
            i = 7;
        C_Misiak_LO[i] = (*(WilsonCoeffsMisiak[LO]))(i);
        C_Misiak_NLO[i] = (*(WilsonCoeffsMisiak[NLO]))(i);
        C_Misiak_NNLO[i] = (*(WilsonCoeffsMisiak[NNLO]))(i);
    }
}

◆ D()

gslpp::complex AmpDB2::D	(	quarks	qq,
		int	k
	)

private

Definition at line 926 of file AmpDB2.cpp.

{
    return cacheD[indexD(qq, k)];
}

◆ delta_1overm()

gslpp::complex AmpDB2::delta_1overm ( quark q )

private

Value of 1/mb corrections of $\Gamma_{21}$ (hep-ph/0612167)

Parameters

[in] quark index of the neutral B mesons @detail requires compute_matrixelements() and Wilson coefficients in Buras basis

Definition at line 985 of file AmpDB2.cpp.

{
    // hep-ph/0612167 equation (10)
    gslpp::complex lambda_c, lambda_u;
    switch (q)
    {
    case d:
        lambda_c = lambda_c_d;
        lambda_u = lambda_u_d;
        break;
    case s:
        lambda_c = lambda_c_s;
        lambda_u = lambda_u_s;
        break;
    default:
        throw std::runtime_error("AmpDB2::delta_1overm(quark q): invalid quark index: ");
    }
    compute_deltas_1overm(q);
    return -(lambda_c * lambda_c * deltas_1overm(cc, q) + 2. * lambda_c * lambda_u * deltas_1overm(cu, q) + lambda_u * lambda_u * deltas_1overm(uu, q)).conjugate();
}

◆ delta_1overm_tradBasis()

gslpp::complex AmpDB2::delta_1overm_tradBasis ( quark q )

private

Value of 1/mb corrections of $\Gamma_{21}$ (hep-ph/0308029v2)

Parameters

[in] quark index of the neutral B mesons @detail requires computeCKMandMasses() before use

Definition at line 945 of file AmpDB2.cpp.

{
    // hep-ph/0308029: equation 18
    compute_deltas_1overm_tradBasis(q);
    switch (q)
    {
    case d:
        return VtbVtd2 * deltas_1overm_tradBasis(uu, d) + 2. * VcbVcd * VtbVtd * (deltas_1overm_tradBasis(uu, d) - deltas_1overm_tradBasis(cu, d)) + VcbVcd2 * (deltas_1overm_tradBasis(cc, d) + deltas_1overm_tradBasis(uu, d) - 2. * deltas_1overm_tradBasis(cu, d));
    case s:
        return VtbVts2 * deltas_1overm_tradBasis(uu, s) + 2. * VcbVcs * VtbVts * (deltas_1overm_tradBasis(uu, s) - deltas_1overm_tradBasis(cu, s)) + VcbVcs2 * (deltas_1overm_tradBasis(cc, s) + deltas_1overm_tradBasis(uu, s) - 2. * deltas_1overm_tradBasis(cu, s));
    default:
        throw std::runtime_error("AmpDB2::delta_1overm_tradBasis(quark q): invalid quark index: ");
    }
}

◆ deltas_1overm()

gslpp::complex AmpDB2::deltas_1overm	(	quarks	qq,
		quark	q
	)

private

Definition at line 974 of file AmpDB2.cpp.

{
    return cache_deltas_1overm[index_deltas(qq, q)];
}

◆ deltas_1overm_tradBasis()

gslpp::complex AmpDB2::deltas_1overm_tradBasis	(	quarks	qq,
		quark	q
	)

private

Definition at line 969 of file AmpDB2.cpp.

{
    return cache_deltas_1overm_tradBasis[index_deltas(qq, q)];
}

◆ F()

double AmpDB2::F	(	quarks	qq,
		int	k,
		int	i,
		int	j
	)

private

Definition at line 753 of file AmpDB2.cpp.

{
    return F0(qq, k, i, j) + as_4pi_mu1 * F1(qq, k, i, j);
}

◆ F0()

double AmpDB2::F0	(	quarks	qq,
		int	k,
		int	i,
		int	j
	)

private

Definition at line 743 of file AmpDB2.cpp.

{
    return cacheF0[indexF(qq, k, i, j)];
}

◆ F1()

double AmpDB2::F1	(	quarks	qq,
		int	k,
		int	i,
		int	j
	)

private

Definition at line 748 of file AmpDB2.cpp.

{
    return cacheF1[indexF(qq, k, i, j)];
}

◆ g()

gslpp::complex AmpDB2::g	(	quarks	qq,
		int	i
	)

private

Definition at line 1063 of file AmpDB2.cpp.

{
    return cacheg[indexg(qq, i)];
}

◆ Gamma21overM21()

gslpp::complex AmpDB2::Gamma21overM21	(	orders	order,
		mass_schemes	mass_scheme,
		int	BMeson
	)

protected

A method to compute $\frac{\Gamma_{21},M_{21}}^{bq}$ @detail source: Marvin Gerlach (2205.07907 and thesis) with 1/mb corrections from Lenz (hep-ph/0612167)

Parameters

[in]	order	the QCD order of the computation
[in]	mass_scheme	the scheme for the bottom quark mass
[in]	BMeson	the B meson index (0 for Bd and 1 for Bs)

Returns: $\frac{\Gamma_{21},M_{21}}^{bq}$

Definition at line 1082 of file AmpDB2.cpp.

{
    computeCKMandMasses(NNLO, mass_scheme);
 
    // calculate M_21_q / <O_1>
    gslpp::complex M21overme0;
    double MBq;
    quark q;
    switch (BMeson)
    {
    case 0:
        q = d;
        // calculate DB=2 matrix elements for usage of "me" and "delta_1overm_(quark)"
        compute_matrixelements(q, FULLNLO);
        M21overme0 = M21_Bd(FULLNLO) * HCUT / me(0);
        MBq = MB;
        break;
    case 1:
        q = s;
        // calculate DB=2 matrix elements for usage of "me" and "delta_1overm_(quark)"
        compute_matrixelements(q, FULLNLO);
        M21overme0 = M21_Bs(FULLNLO) * HCUT / me(0);
        MBq = MB_s;
        break;
    default:
        throw std::runtime_error("AmpDB2::Gamma21overM21(orders order, int BMeson): invalid B meson index: ");
    }
 
    // calculate DB=1 Wilson coefficients
    computeWilsonCoeffsMisiak();
 
    // calculate DB=2 coefficients in pole scheme for usage of "c_H()"
    compute_pp_s();
 
    // switch to another scheme if needed
    if (mass_scheme == MSbar)
    {
        poletoMSbar_pp_s();
    }
    else if (mass_scheme == PS)
    {
        poletoPS_pp_s();
    }
    else if (mass_scheme == pole or mass_scheme == only1overmb)
    {
    }
    else if (mass_scheme == MSbar_partialNNLO)
    {
        compute_partialNNLO();
        poletoMSbar_pp_s();
    }
    else if (mass_scheme == PS_partialNNLO)
    {
        compute_partialNNLO();
        poletoPS_pp_s();
    }
    else if (mass_scheme == MSbar_partialN3LO)
    {
        poletoMSbar_pp_s_partialN3LO();
    }
    else if (mass_scheme == PS_partialN3LO)
    {
        poletoPS_pp_s_partialN3LO();
    }
    else
    {
        std::cerr << "WARNING: mass_scheme might no be implemented.\n";
    }
 
    /* partial contribution without 1overm
    gslpp::vector<gslpp::complex> Gamma21overM21_partial = Mb2_prefactor * (c_H_partial(q, 0) + c_H_partial(q, 1) * me(1)/me(0) + c_H_partial(q, 2) * me(2)/me(0)).conjugate()
            * Gf2 / (24 * M_PI * MBq) / M21overme0; */
 
    // Gerlach thesis eq. 6.1 divided by M_21
    gslpp::complex Gamma21overM21 = 0.;
    if (flag_RI)
    {
        Gamma21overM21 = Mb2_prefactor * (c_H(q, LO) * meoverme0).conjugate();
        // compute_matrixelements(q, LO); //minimal difference
        Gamma21overM21 += Mb2_prefactor * (c_H(q, NLO) * meoverme0).conjugate();
    }
    else
    {
        Gamma21overM21 = Mb2_prefactor * (c_H(q, order) * meoverme0).conjugate();
    }
    computeWilsonCoeffsBuras(); /*calculate DB=1 Wilson coefficients in the basis for "delta_1overm" */
    Gamma21overM21 += Mb2_prefactor_1overm * delta_1overm(q) / me(0);
    Gamma21overM21 *= Gf2 / (24. * M_PI * MBq) / M21overme0;
    return Gamma21overM21;
}

◆ Gamma21overM21_tradBasis()

gslpp::complex AmpDB2::Gamma21overM21_tradBasis	(	orders	order,
		quark	q
	)

protected

A method to compute $\frac{\Gamma_{21},M_{21}}^{bq}$ in the traditional basis @detail source: Ciuchini (hep-ph/0308029v2)

Parameters

[in]	order	the QCD order of the computation
[in]	BMeson	the B meson index (0 for Bd and 1 for Bs)

Returns: $\frac{\Gamma_{21},M_{21}}^{bq}$

Definition at line 256 of file AmpDB2.cpp.

{
    computeCKMandMasses(NLO);
   // calculate DB=2 matrix elements for usage of "me" and "delta_1overm_tradBasis(quark)"
    compute_matrixelements(q, FULLNLO);
 
    // calculate M_21_q / <O_1>
    gslpp::complex M21overme0;
    double MBq;
    switch (q)
    {
    case d:
        M21overme0 = M21_Bd(FULLNLO) * HCUT / me(0); 
        MBq = MB;
        break;
    case s:
        M21overme0 = M21_Bs(FULLNLO) * HCUT / me(0); 
        MBq = MB_s;
        break;
    default:
        throw std::runtime_error("AmpDB2::Gamma21overM21_tradBasis(orders order, quark q): invalid quark index: ");
    }
    
    // calculate DB=1 Wilson coefficients
    computeWilsonCoeffsBuras();
 
    // calculate DB=2 coefficients for usage of "c(quark)"
    computeF0();
    computeF1();
    computeP();
 
    // staying in the old basis: hep-ph/0308029v2: eq. 16 divided by M_21
    /*
    computeD(order);
    gslpp::complex Gamma21overM21 = -Gf2 / (24 * M_PI * MBq) / M21overme0 *
                (Mb2_prefactor * (c(q)(0) + c(q)(1) * me(1)/me(0) + c(q)(2) * me(2)/me(0)) +
            Mb2_prefactor_1overm * delta_1overm_tradBasis(q)/me(0));
    */
 
    // switching to the new basis
    gslpp::complex Gamma21overM21;
    switch (order)
    {
    case FULLNLO:
        computeD(FULLNLO);
        Gamma21overM21 = Mb2_prefactor * (c(q, LO) * meoverme0) + Mb2_prefactor_1overm * (-delta_1overm(q)) / me(0);
        computeD(LO);
        Gamma21overM21 += Mb2_prefactor * (c(q, NLO) * meoverme0);
        Gamma21overM21 *= -Gf2 / (24. * M_PI * MBq) / M21overme0;
        break;
    case LO:
        computeD(LO);
        Gamma21overM21 = -Gf2 / (24. * M_PI * MBq) / M21overme0 *
                         (Mb2_prefactor * (c(q, LO) * meoverme0) + Mb2_prefactor_1overm * (-delta_1overm(q)) / me(0));
        break;
    default:
        throw std::runtime_error("AmpDB2::Gamma21overM21_tradBasis(orders order, quark q): order not implemented");
    }
    return Gamma21overM21;
}

◆ getGamma21overM21()

gslpp::complex AmpDB2::getGamma21overM21	(	orders	order,
		mass_schemes	mass_scheme = `MSbar`
	)

inline

The value of $\frac{\Gamma_{21},M_{21}}$ from Gerlach (2205.07907 and thesis)

Parameters

[in]	order	the QCD order of the computation
[in]	mass_scheme	the scheme for the bottom quark mass

Returns: $\frac{\Gamma_{21},M_{21}}$

Definition at line 91 of file AmpDB2.h.

    {
        return Gamma21overM21(order, mass_scheme, BMeson);
    }

◆ getGamma21overM21_tradBasis()

gslpp::complex AmpDB2::getGamma21overM21_tradBasis ( orders order )

inline

The value of $\frac{\Gamma_{21},M_{21}}$ in the traditional basis.

Parameters

[in]	order	the QCD order of the computation
[in]	mass_scheme	the scheme for the bottom quark mass

Returns: $\frac{\Gamma_{21},M_{21}}$

Definition at line 73 of file AmpDB2.h.

    {
        if (BMeson == 0)
        {
            return Gamma21overM21_tradBasis(order, d);
        }
        else
        {
            return Gamma21overM21_tradBasis(order, s);
        }
    }

◆ getM21()

gslpp::complex AmpDB2::getM21 ( orders order )

inline

The value of $M_{21}$ for $B_{d,s}$ mesons.

Parameters

[in] order the QCD order of the computation

Returns: $M_{21}$

Definition at line 51 of file AmpDB2.h.

    {
        if (BMeson == 0)
        {
            return M21_Bd(order);
        }
        else
        {
            return M21_Bs(order);
        }
    }

◆ getPB()

gslpp::complex AmpDB2::getPB ( )

inline

Definition at line 96 of file AmpDB2.h.

    {
        if (BMeson == 0)
        {
            return PBd();
        }
        else
        {
            return PBs();
        }
    }

◆ getRB()

gslpp::complex AmpDB2::getRB ( orders order )

inline

Definition at line 107 of file AmpDB2.h.

    {
        if (BMeson == 0)
        {
            return RBd(order);
        }
        else
        {
            return RBs(order);
        }
    }

◆ gtilde()

gslpp::complex AmpDB2::gtilde	(	quarks	qq,
		int	i
	)

private

Definition at line 1068 of file AmpDB2.cpp.

{
    return cachegtilde[indexg(qq, i)];
}

◆ H()

gslpp::complex AmpDB2::H	(	quarks	qq,
		orders	order
	)

private

Definition at line 1232 of file AmpDB2.cpp.

{
    if (order == LO)
        return H_partial(qq, 1, 8, 1, 8, 0);
    if (order == NLO)
        return H_partial(qq, 1, 8, 1, 8, 1);
    if (order == NNLO)
        return H_partial(qq, 1, 8, 1, 8, 2);
    if (order == FULLNLO)
        return H_partial(qq, 1, 8, 1, 8, 0) + H_partial(qq, 1, 8, 1, 8, 1);
    if (order == FULLNNLO)
        return H_partial(qq, 1, 8, 1, 8, 0) + H_partial(qq, 1, 8, 1, 8, 1) + H_partial(qq, 1, 8, 1, 8, 2);
    if (order == FULLNNNLO)
        return H_partial(qq, 1, 8, 1, 8, 0) + H_partial(qq, 1, 8, 1, 8, 1) + H_partial(qq, 1, 8, 1, 8, 2) + H_partial(qq, 1, 8, 1, 8, 3);
    throw std::runtime_error("AmpDB2::H(quarks qq, orders order) order not implemented");
}

◆ H_allpartial()

gslpp::vector< gslpp::complex > AmpDB2::H_allpartial ( quarks qq )

private

Definition at line 1295 of file AmpDB2.cpp.

{
    gslpp::vector<gslpp::complex> result(12, 0.);
    result.assign(0, H_partial(qq, 1, 2, 1, 2, 0));
    result.assign(1, H_partial(qq, 1, 2, 1, 2, 1));
    result.assign(2, H_partial(qq, 1, 2, 1, 2, 2));
    result.assign(3, H_partial(qq, 1, 2, 3, 6, 0));
    result.assign(4, H_partial(qq, 1, 2, 3, 6, 1));
    result.assign(5, H_partial(qq, 3, 6, 3, 6, 0));
    result.assign(6, H_partial(qq, 3, 6, 3, 6, 1));
    result.assign(7, H_partial(qq, 1, 2, 8, 8, 1));
    result.assign(8, H_partial(qq, 1, 2, 8, 8, 2));
    result.assign(9, H_partial(qq, 3, 6, 8, 8, 1));
    result.assign(10, H_partial(qq, 3, 6, 8, 8, 2));
    result.assign(11, H_partial(qq, 8, 8, 8, 8, 2));
    return result;
}

◆ H_partial()

gslpp::complex AmpDB2::H_partial	(	quarks	qq,
		int	i_start,
		int	i_end,
		int	j_start,
		int	j_end,
		int	n
	)

private

Definition at line 1331 of file AmpDB2.cpp.

{
    gslpp::complex result = 0.;
    for (int i = i_start; i <= i_end; i++)
    {
        if (i == 7)
            i++;
        for (int j = j_start; j <= j_end; j++)
        {
            if (j == 7)
                j++;
            if (j < i)
                continue;
            if (n == 0)
            {
                result += C_Misiak_LO[i - 1] * C_Misiak_LO[j - 1] * p(qq, i, j, 0);
            }
            else if (n == 1)
            {
                if (j == 1 or j == 2 or j == 8)
                {
                    result += as_4pi_mu1 * C_Misiak_LO[i - 1] * C_Misiak_LO[j - 1] * p(qq, i, j, 1) + (C_Misiak_NLO[i - 1] * C_Misiak_LO[j - 1] + C_Misiak_LO[i - 1] * C_Misiak_NLO[j - 1]) * p(qq, i, j, 0);
                }
                else if (3 <= j and j <= 6)
                {
                    result += as_4pi_mu1 * C_Misiak_LO[i - 1] * C_Misiak_LO[j - 1] * p(qq, i, j, 1) + (C_Misiak_NLO[i - 1] * C_Misiak_LO[j - 1] + C_Misiak_LO[i - 1] * C_Misiak_NLO[j - 1]) * p(qq, i, j, 0, flag_LOz);
                }
            }
            else if (n == 2)
            {
                if (j == 1 or j == 2 or j == 8)
                {
                    result += as_4pi_mu1 * as_4pi_mu1 * C_Misiak_LO[i - 1] * C_Misiak_LO[j - 1] * p(qq, i, j, 2) + as_4pi_mu1 * (C_Misiak_NLO[i - 1] * C_Misiak_LO[j - 1] + C_Misiak_LO[i - 1] * C_Misiak_NLO[j - 1]) * p(qq, i, j, 1, flag_LOz) + (C_Misiak_NNLO[i - 1] * C_Misiak_LO[j - 1] + C_Misiak_NLO[i - 1] * C_Misiak_NLO[j - 1] + C_Misiak_LO[i - 1] * C_Misiak_NNLO[j - 1]) * p(qq, i, j, 0, flag_LOz);
                }
            }
            else if (n == 3)
            {
                result += as_4pi_mu1 * as_4pi_mu1 * as_4pi_mu1 * C_Misiak_LO[i - 1] * C_Misiak_LO[j - 1] * p(qq, i, j, 3) + as_4pi_mu1 * as_4pi_mu1 * (C_Misiak_NLO[i - 1] * C_Misiak_LO[j - 1] + C_Misiak_LO[i - 1] * C_Misiak_NLO[j - 1]) * p(qq, i, j, 2) + as_4pi_mu1 * (C_Misiak_NNLO[i - 1] * C_Misiak_LO[j - 1] + C_Misiak_NLO[i - 1] * C_Misiak_NLO[j - 1] + C_Misiak_LO[i - 1] * C_Misiak_NNLO[j - 1]) * p(qq, i, j, 1) + (C_Misiak_NNLO[i - 1] * C_Misiak_NLO[j - 1] + C_Misiak_NLO[i - 1] * C_Misiak_NNLO[j - 1]) * p(qq, i, j, 0);
            }
            else
            {
                throw(std::runtime_error("AmpDB2::H_partial order not implemented"));
            }
        }
    }
    return result;
}

◆ H_s()

gslpp::complex AmpDB2::H_s	(	quarks	qq,
		orders	order
	)

private

Definition at line 1248 of file AmpDB2.cpp.

{
    if (order == LO)
        return H_s_partial(qq, 1, 8, 1, 8, 0);
    if (order == NLO)
        return H_s_partial(qq, 1, 8, 1, 8, 1);
    if (order == NNLO)
        return H_s_partial(qq, 1, 8, 1, 8, 2);
    if (order == FULLNLO)
        return H_s_partial(qq, 1, 8, 1, 8, 0) + H_s_partial(qq, 1, 8, 1, 8, 1);
    if (order == FULLNNLO)
        return H_s_partial(qq, 1, 8, 1, 8, 0) + H_s_partial(qq, 1, 8, 1, 8, 1) + H_s_partial(qq, 1, 8, 1, 8, 2);
    if (order == FULLNNNLO)
        return H_s_partial(qq, 1, 8, 1, 8, 0) + H_s_partial(qq, 1, 8, 1, 8, 1) + H_s_partial(qq, 1, 8, 1, 8, 2) + H_s_partial(qq, 1, 8, 1, 8, 3);
    throw std::runtime_error("AmpDB2::H_s(quarks qq, orders order) order not implemented");
}

◆ H_s_allpartial()

gslpp::vector< gslpp::complex > AmpDB2::H_s_allpartial ( quarks qq )

private

Definition at line 1313 of file AmpDB2.cpp.

{
    gslpp::vector<gslpp::complex> result(12, 0.);
    result.assign(0, H_s_partial(qq, 1, 2, 1, 2, 0));
    result.assign(1, H_s_partial(qq, 1, 2, 1, 2, 1));
    result.assign(2, H_s_partial(qq, 1, 2, 1, 2, 2));
    result.assign(3, H_s_partial(qq, 1, 2, 3, 6, 0));
    result.assign(4, H_s_partial(qq, 1, 2, 3, 6, 1));
    result.assign(5, H_s_partial(qq, 3, 6, 3, 6, 0));
    result.assign(6, H_s_partial(qq, 3, 6, 3, 6, 1));
    result.assign(7, H_s_partial(qq, 1, 2, 8, 8, 1));
    result.assign(8, H_s_partial(qq, 1, 2, 8, 8, 2));
    result.assign(9, H_s_partial(qq, 3, 6, 8, 8, 1));
    result.assign(10, H_s_partial(qq, 3, 6, 8, 8, 2));
    result.assign(11, H_s_partial(qq, 8, 8, 8, 8, 2));
    return result;
}

◆ H_s_partial()

gslpp::complex AmpDB2::H_s_partial	(	quarks	qq,
		int	i_start,
		int	i_end,
		int	j_start,
		int	j_end,
		int	n
	)

private

Definition at line 1379 of file AmpDB2.cpp.

{
    gslpp::complex result = 0.;
    for (int i = i_start; i <= i_end; i++)
    {
        if (i == 7)
            i++;
        for (int j = j_start; j <= j_end; j++)
        {
            if (j == 7)
                j++;
            if (j < i)
                continue;
            if (n == 0)
            {
                result += C_Misiak_LO[i - 1] * C_Misiak_LO[j - 1] * p_s(qq, i, j, 0);
            }
            else if (n == 1)
            {
                if (j == 1 or j == 2 or j == 8)
                {
                    result += as_4pi_mu1 * C_Misiak_LO[i - 1] * C_Misiak_LO[j - 1] * p_s(qq, i, j, 1) + (C_Misiak_NLO[i - 1] * C_Misiak_LO[j - 1] + C_Misiak_LO[i - 1] * C_Misiak_NLO[j - 1]) * p_s(qq, i, j, 0);
                }
                else if (3 <= j and j <= 6)
                {
                    result += as_4pi_mu1 * C_Misiak_LO[i - 1] * C_Misiak_LO[j - 1] * p_s(qq, i, j, 1) + (C_Misiak_NLO[i - 1] * C_Misiak_LO[j - 1] + C_Misiak_LO[i - 1] * C_Misiak_NLO[j - 1]) * p_s(qq, i, j, 0, flag_LOz);
                }
            }
            else if (n == 2)
            {
                if (j == 1 or j == 2 or j == 8)
                {
                    result += as_4pi_mu1 * as_4pi_mu1 * C_Misiak_LO[i - 1] * C_Misiak_LO[j - 1] * p_s(qq, i, j, 2) + as_4pi_mu1 * (C_Misiak_NLO[i - 1] * C_Misiak_LO[j - 1] + C_Misiak_LO[i - 1] * C_Misiak_NLO[j - 1]) * p_s(qq, i, j, 1, flag_LOz) + (C_Misiak_NNLO[i - 1] * C_Misiak_LO[j - 1] + C_Misiak_NLO[i - 1] * C_Misiak_NLO[j - 1] + C_Misiak_LO[i - 1] * C_Misiak_NNLO[j - 1]) * p_s(qq, i, j, 0, flag_LOz);
                }
            }
            else if (n == 3)
            {
                result += as_4pi_mu1 * as_4pi_mu1 * as_4pi_mu1 * C_Misiak_LO[i - 1] * C_Misiak_LO[j - 1] * p_s(qq, i, j, 3) + as_4pi_mu1 * as_4pi_mu1 * (C_Misiak_NLO[i - 1] * C_Misiak_LO[j - 1] + C_Misiak_LO[i - 1] * C_Misiak_NLO[j - 1]) * p_s(qq, i, j, 2) + as_4pi_mu1 * (C_Misiak_NNLO[i - 1] * C_Misiak_LO[j - 1] + C_Misiak_NLO[i - 1] * C_Misiak_NLO[j - 1] + C_Misiak_LO[i - 1] * C_Misiak_NNLO[j - 1]) * p_s(qq, i, j, 1) + (C_Misiak_NNLO[i - 1] * C_Misiak_NLO[j - 1] + C_Misiak_NLO[i - 1] * C_Misiak_NNLO[j - 1]) * p_s(qq, i, j, 0);
            }
            else
            {
                throw(std::runtime_error("AmpDB2::H_s_partial order not implemented"));
            }
        }
    }
    return result;
}

◆ index_deltas()

int AmpDB2::index_deltas	(	quarks	qq,
		quark	q
	)

private

Definition at line 980 of file AmpDB2.cpp.

{
    return qq * 2 + q;
}

◆ index_p()

int AmpDB2::index_p	(	quarks	qq,
		int	i,
		int	j,
		int	n
	)

private

Definition at line 1452 of file AmpDB2.cpp.

{
    return n * 192 + qq * 64 + (i - 1) * 8 + (j - 1);
}

◆ indexD()

int AmpDB2::indexD	(	quarks	qq,
		int	k
	)

private

Definition at line 932 of file AmpDB2.cpp.

{
    if (k != 1 and k != 2)
    {
        throw std::runtime_error("AmpDB2::indexD(quarks qq, int k): invalid k");
    }
    return qq * 2 + (k - 1);
}

◆ indexF()

int AmpDB2::indexF	(	quarks	qq,
		int	k,
		int	i,
		int	j
	)

private

Definition at line 764 of file AmpDB2.cpp.

{
    return qq * 8 + (k - 1) * 4 + (i - 1) * 2 + (j - 1);
}

◆ indexg()

int AmpDB2::indexg	(	quarks	qq,
		int	i
	)

private

Definition at line 1073 of file AmpDB2.cpp.

{
    return qq * 4 + i;
}

◆ indexP()

int AmpDB2::indexP	(	quarks	qq,
		int	k,
		int	i,
		int	j
	)

private

Definition at line 769 of file AmpDB2.cpp.

{
    return qq * 28 + (k - 1) * 14 + (i - 1) * 7 + (j - 2);
}

◆ M21_Bd()

gslpp::complex AmpDB2::M21_Bd ( orders order )

protected

A method to compute $M_{21}^{bd}$.

Parameters

[in] order the QCD order of the computation

Returns: $M_{21}^{bd}$

Definition at line 163 of file AmpDB2.cpp.

{
    if (mySM.getFlavour().getHDF2().getCoeffBd().getOrder() < getHighest(order))
        throw std::runtime_error("DmBd::computeThValue(): requires cofficient of order not computed");
 
    gslpp::vector<gslpp::complex> **allcoeff = mySM.getFlavour().ComputeCoeffBd(
        mySM.getBBd().getMu(),
        mySM.getBBd().getScheme());
 
    gslpp::vector<double> me(mySM.getBBd().getBpars());
    double MBd = mySM.getMesons(QCD::B_D).getMass();
    double Mb = mySM.Mrun(mySM.getBBd().getMu(),
                          mySM.getQuarks(QCD::BOTTOM).getMass_scale(),
                          mySM.getQuarks(QCD::BOTTOM).getMass(), QCD::BOTTOM, FULLNNLO);
    double Md = mySM.Mrun(mySM.getBBd().getMu(),
                          mySM.getQuarks(QCD::DOWN).getMass_scale(),
                          mySM.getQuarks(QCD::DOWN).getMass(), QCD::DOWN, FULLNNLO);
    double KBd = MBd / (Mb + Md) * MBd / (Mb + Md);
    double Fb = mySM.getMesons(QCD::B_D).getDecayconst();
    me(0) *= 1. / 3. * MBd * Fb * Fb;
    me(1) *= -5. / 24. * KBd * MBd * Fb * Fb;
    me(2) *= 1. / 24. * KBd * MBd * Fb * Fb;
    me(3) *= 1. / 4. * KBd * MBd * Fb * Fb;
    me(4) *= 1. / 12. * KBd * MBd * Fb * Fb;
 
#if SUSYFIT_DEBUG & 1
    std::cout << "Bd: me(0) = " << me(0) << std::endl;
#endif
#if SUSYFIT_DEBUG & 2
    std::cout << "coefficient Bd: " << (*(allcoeff[LO]) + *(allcoeff[NLO]))(0) << std::endl;
    std::cout << "M: " << me << std::endl;
    std::cout << "mu : " << mySM.getBBd().getMu() << ", mut: " << mySM.getMut() << ", scheme: " << mySM.getBBd().getScheme() << ", B par.: " << mySM.getBBd().getBpars()(0) << std::endl;
    std::cout << "U (mut): " << (mySM.getFlavour().getHDF2().getUDF2().Df2Evol(mySM.getBBd().getMu(), mySM.getMut(), LO)(0, 0) + mySM.getFlavour().getHDF2().getUDF2().Df2Evol(mySM.getBBd().getMu(), mySM.getMut(), NLO)(0, 0)) << std::endl;
#endif
 
    switch (order)
    {
    case FULLNLO:
        return ((*(allcoeff[LO]) + *(allcoeff[NLO])) * me / HCUT);
    case LO:
        return ((*(allcoeff[LO])) * me / HCUT);
    default:
        throw std::runtime_error("AmpDB2::AmpBd(): order not implemented");
    }
}

◆ M21_Bs()

gslpp::complex AmpDB2::M21_Bs ( orders order )

protected

A method to compute $M_{21}^{bs}$.

Parameters

[in] order the QCD order of the computation

Returns: $M_{21}^{bs}$

Definition at line 209 of file AmpDB2.cpp.

{
    if (mySM.getFlavour().getHDF2().getCoeffBs().getOrder() < getHighest(order))
        throw std::runtime_error("DmBd::computeThValue(): requires cofficient of order not computed");
 
    // Wilson coefficients in same mass scale and scheme as B parameters
    gslpp::vector<gslpp::complex> **allcoeff = mySM.getFlavour().ComputeCoeffBs(
        mySM.getBBs().getMu(),
        mySM.getBBs().getScheme());
 
    gslpp::vector<double> me(mySM.getBBs().getBpars());
    double MBs = mySM.getMesons(QCD::B_S).getMass();
    double Mb = mySM.Mrun(mySM.getBBs().getMu(),
                          mySM.getQuarks(QCD::BOTTOM).getMass_scale(),
                          mySM.getQuarks(QCD::BOTTOM).getMass(), QCD::BOTTOM, FULLNNLO);
    double Ms = mySM.Mrun(mySM.getBBs().getMu(),
                          mySM.getQuarks(QCD::STRANGE).getMass_scale(),
                          mySM.getQuarks(QCD::STRANGE).getMass(), QCD::STRANGE, FULLNNLO);
    double KBs = MBs / (Mb + Ms) * MBs / (Mb + Ms);
    double Fbs = mySM.getMesons(QCD::B_S).getDecayconst();
 
    // matrix elements as in hep-ph/9604387v2 with KBs independent terms in 4 and 5
    // differ by factor of 2 to arXiv:1602.03560v2 etc.
    me(0) *= 1. / 3. * MBs * Fbs * Fbs;
    me(1) *= -5. / 24. * KBs * MBs * Fbs * Fbs;
    me(2) *= 1. / 24. * KBs * MBs * Fbs * Fbs;
    me(3) *= 1. / 4. * KBs * MBs * Fbs * Fbs;
    me(4) *= 1. / 12. * KBs * MBs * Fbs * Fbs;
#if SUSYFIT_DEBUG & 1
    std::cout << "Bs: me(0) = " << me(0) << std::endl;
#endif
 
    switch (order)
    {
    case FULLNLO:
        return ((*(allcoeff[LO]) + *(allcoeff[NLO])) * me / HCUT);
    case LO:
        return ((*(allcoeff[LO])) * me / HCUT);
    default:
        throw std::runtime_error("AmpDB2::AmpBs(): order not implemented");
    }
}

◆ p()

double AmpDB2::p	(	quarks	qq,
		int	i,
		int	j,
		int	n,
		bool	flag_LOz = `false`
	)

private

Definition at line 1427 of file AmpDB2.cpp.

{
    if (n < 0 or n >= 768)
        throw(std::runtime_error("AmpDB2::p(quarks qq, int i, int j, int n, bool flag_LOz) out of index"));
    if (flag_LOz)
    {
        if (n >= 576)
            throw(std::runtime_error("AmpDB2::p(quarks qq, int i, int j, int n, bool flag_LOz) out of index"));
        return cache_p_LO[index_p(qq, i, j, n)];
    }
    return cache_p[index_p(qq, i, j, n)];
}

◆ P()

double AmpDB2::P	(	quarks	qq,
		int	k,
		int	i,
		int	j
	)

private

Definition at line 758 of file AmpDB2.cpp.

{
    return cacheP[indexP(qq, k, i, j)];
}

◆ p_s()

double AmpDB2::p_s	(	quarks	qq,
		int	i,
		int	j,
		int	n,
		bool	flag_LOz = `false`
	)

private

Definition at line 1439 of file AmpDB2.cpp.

{
    if (n < 0 or n >= 768)
        throw(std::runtime_error("AmpDB2::p_s(quarks qq, int i, int j, int n, bool flag_LOz) out of index"));
    if (flag_LOz)
    {
        if (n >= 576)
            throw(std::runtime_error("AmpDB2::p_s(quarks qq, int i, int j, int n, bool flag_LOz) out of index"));
        return cache_ps_LO[index_p(qq, i, j, n)];
    }
    return cache_ps[index_p(qq, i, j, n)];
}

◆ PBd()

gslpp::complex AmpDB2::PBd ( )

protected

Definition at line 2589 of file AmpDB2.cpp.

{
    double mbpole = mySM.Mbar2Mp(mySM.getQuarks(QCD::BOTTOM).getMass(), QCD::BOTTOM);
    double Mw = mySM.Mw();
    double kappa = -2. * M_PI * mbpole * mbpole /
                   (3. * Mw * Mw * mySM.getFlavour().getHDF2().getUDF2().etabS0(mySM.getBBd().getMu()));
 
    double n[13] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
 
    n[0] = 0.1797;
    n[1] = 0.1391;
    n[5] = 1.0116;
    n[6] = 0.0455;
    n[8] = -0.0714;
    n[10] = -0.3331;
 
    double B1 = mySM.getBBd().getBpars()(0);
    double B2 = mySM.getBBd().getBpars()(1);
 
    gslpp::complex PBd = -2. * kappa / mySM.getCBd() *
                         (gslpp::complex(1, 2. * mySM.getPhiBd(), true) * (n[0] + (n[5] * B2 + n[10]) / B1) - gslpp::complex(1. / mySM.getCKM().computeRt(), mySM.getCKM().computeBeta() + 2. * mySM.getPhiBd(), true) * (n[1] + (n[6] * B2 + n[11]) / B1) + gslpp::complex(1. / mySM.getCKM().computeRt() / mySM.getCKM().computeRt(), 2. * (mySM.getCKM().computeBeta() + mySM.getPhiBd()), true) * (n[2] + (n[7] * B2 + n[12]) / B1));
 
    return PBd;
}

◆ PBs()

gslpp::complex AmpDB2::PBs ( )

protected

Definition at line 2614 of file AmpDB2.cpp.

{
    double mbpole = mySM.Mbar2Mp(mySM.getQuarks(QCD::BOTTOM).getMass(), QCD::BOTTOM);
 
    double Mw = mySM.Mw();
    double kappa = -2. * M_PI * mbpole * mbpole /
                   (3. * Mw * Mw * mySM.getFlavour().getHDF2().getUDF2().etabS0(mySM.getBBs().getMu()));
 
    double n[13] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
 
    n[0] = 0.1797;
    n[1] = 0.1391;
    n[5] = 1.0116;
    n[6] = 0.0455;
    n[8] = -0.0714;
    n[10] = -0.3331;
 
    double B1 = mySM.getBBs().getBpars()(0);
    double B2 = mySM.getBBs().getBpars()(1);
 
    gslpp::complex PBs = -2. * kappa / mySM.getCBs() *
                         (gslpp::complex(1, 2. * mySM.getPhiBs(), true) * (n[0] + (n[5] * B2 + n[10]) / B1) - gslpp::complex(1. / mySM.getCKM().computeRts(), -mySM.getCKM().computeBetas() + 2. * mySM.getPhiBs(), true) * (n[1] + (n[6] * B2 + n[11]) / B1) + gslpp::complex(1. / mySM.getCKM().computeRts() / mySM.getCKM().computeRts(), 2. * (-mySM.getCKM().computeBetas() + mySM.getPhiBs()), true) * (n[2] + (n[7] * B2 + n[12]) / B1));
 
    return PBs;
}

◆ poletoMSbar_pp_s()

void AmpDB2::poletoMSbar_pp_s ( )

private

Definition at line 2324 of file AmpDB2.cpp.

{
    // arxiv:2106.05979 eq. (33)
    double log_mub_Mb = 2. * log(mu_b / Mb);
    PoletoMS_as1 = 4. / 3. + log_mub_Mb;
    PoletoMS_as2_z0 = (307. / 32. + M_PI2 / 3. + M_PI2 / 9. * log2 - 1. / 6. * zeta3 - (71. / 144. + M_PI2 / 18.) * nl + log_mub_Mb * (493. / 72. - nl * 13. / 36. + log_mub_Mb * (43. / 24. - nl / 12.)));
    PoletoMS_as2_z1 = 4. / 3. * M_PI2 / 8. * Mc / Mb - z;
    // a_s(mu_b) to a_s(mu_1)
    double beta0 = 23. / 3.;
    double log_mu1_mub = 2. * log(mu_1 / mu_b);
    double as_running1 = beta0 * log_mu1_mub;
    bool flag_LOz = true;
    for (quarks qq = cc; qq <= uu; qq = quarks(qq + 1))
    {
        for (int i = 1; i <= 6; i++)
        {
            // not all terms used for n=2
            for (int j = i; j <= 8; j++)
            {
                if (j == 3)
                    j = 8;
                if (i >= 3)
                    j = 8;
                cache_p[index_p(qq, i, j, 2)] += 8. * PoletoMS_as1 * p(qq, i, j, 1, flag_LOz) + (32. * PoletoMS_as2_z0 + 16. * PoletoMS_as1 * PoletoMS_as1 + 8. * PoletoMS_as1 * as_running1) * p(qq, i, j, 0, flag_LOz) + 16. * 2. * PoletoMS_as2_z1 * p(uu, i, j, 0, flag_LOz);
                cache_ps[index_p(qq, i, j, 2)] += 8. * PoletoMS_as1 * p_s(qq, i, j, 1, flag_LOz) + (32. * PoletoMS_as2_z0 + 16. * PoletoMS_as1 * PoletoMS_as1 + 8. * PoletoMS_as1 * as_running1) * p_s(qq, i, j, 0, flag_LOz) + 16. * 2. * PoletoMS_as2_z1 * p_s(uu, i, j, 0, flag_LOz);
                cache_p_LO[index_p(qq, i, j, 2)] = cache_p[index_p(qq, i, j, 2)];
                cache_ps_LO[index_p(qq, i, j, 2)] = cache_ps[index_p(qq, i, j, 2)];
            }
            for (int j = i; j <= 8; j++)
            {
                cache_p_LO[index_p(qq, i, j, 1)] += 8. * PoletoMS_as1 * p(qq, i, j, 0, flag_LOz);
                cache_ps_LO[index_p(qq, i, j, 1)] += 8. * PoletoMS_as1 * p_s(qq, i, j, 0, flag_LOz);
                if (j == 1 or j == 2 or j == 8)
                {
                    cache_p[index_p(qq, i, j, 1)] += 8. * PoletoMS_as1 * p(qq, i, j, 0);
                    cache_ps[index_p(qq, i, j, 1)] += 8. * PoletoMS_as1 * p_s(qq, i, j, 0);
                }
                else if (3 <= j and j <= 6)
                {
                    cache_p[index_p(qq, i, j, 1)] += 8. * PoletoMS_as1 * p(qq, i, j, 0, flag_LOz);
                    cache_ps[index_p(qq, i, j, 1)] += 8. * PoletoMS_as1 * p_s(qq, i, j, 0, flag_LOz);
                }
            }
        }
    }
    return;
}

◆ poletoMSbar_pp_s_partialN3LO()

void AmpDB2::poletoMSbar_pp_s_partialN3LO ( )

private

Definition at line 2426 of file AmpDB2.cpp.

{
    // arxiv:2106.05979 eq. (33)
    double log_mub_Mb = 2. * log(mu_b / Mb);
    double log_mub_Mb_2 = log_mub_Mb * log_mub_Mb;
    double log_mub_Mb_3 = log_mub_Mb_2 * log_mub_Mb;
    double nl2 = nl * nl;
    double Dilogsqrtz = gslpp_special_functions::dilog(sqrtz);
    double Dilogminussqrtz = gslpp_special_functions::dilog(-sqrtz);
    double oneminussqrtz2 = (1. - sqrtz) * (1. - sqrtz);
    double log1minussqrtz = log(1. - sqrtz);
    double log1plussqrtz = log(1. + sqrtz);
    double sqrtz3 = sqrtz * sqrtz * sqrtz;
    double log_mub_Mc = 2. * log(mu_b / Mc);
    PoletoMS_as1 = 4. / 3. + log_mub_Mb;
    PoletoMS_as2 = (307. / 32. + M_PI2 / 3. + M_PI2 / 9. * log2 - 1. / 6. * zeta3 - (71. / 144. + M_PI2 / 18.) * nl + log_mub_Mb * (493. / 72. - nl * 13. / 36. + log_mub_Mb * (43. / 24. - nl / 12.))) + 4. / 3. * M_PI2 / 8. * sqrtz - z;
    PoletoMS_as3 = (8481925. / 93312. + (2. * Dilogminussqrtz) / 9. + (2. * Dilogsqrtz) / 9. - (55. * log2_4) / 162. + (652841. * M_PI2) / 38880. - (575. * log2 * M_PI2) / 162. - (22. * log2_2 * M_PI2) / 81. - (695. * M_PI4) / 7776. + (137. * nl) / 216. - (M_PI2 * nl) / 27. +
                    (log_mub_Mb_3 * (1591. - 160. * nl + 4. * nl2)) / 432. + (log_mub_Mb_2 * (19315. - 2206. * nl + 52. * nl2)) / 864. - (220. * polylog4_12) / 27. + 22.968067156 * sqrtz + (2. * Dilogminussqrtz * sqrtz) / 9. - (2. * Dilogsqrtz * sqrtz) / 9. - 0.982666666 * nl * sqrtz +
                    1.022111111 * sqrtz3 + (2. * Dilogminussqrtz * sqrtz3) / 9. - (2. * Dilogsqrtz * sqrtz3) / 9. + (M_PI2 * sqrtz3) / 9. + 4.014666667 * z - 0.300333333 * nl * z - (log_mub_Mc * sqrtz * (3. * log1minussqrtz * logz - 3. * log1plussqrtz * logz - 3. * M_PI2 + 24. * sqrtz + 6. * logz * sqrtz + 6. * log2z * sqrtz3 - 12. * log1minussqrtz * logz * sqrtz3 - 12. * log1plussqrtz * logz * sqrtz3 + 4. * M_PI2 * sqrtz3 + 9. * log1minussqrtz * logz * z - 9. * log1plussqrtz * logz * z - 9. * M_PI2 * z - 6. * Dilogminussqrtz * (1. + 4. * sqrtz3 + 3. * z) + Dilogsqrtz * (6. - 24. * sqrtz3 + 18. * z))) / 18. + 0.0493333333 * zeta2 + (2. * Dilogminussqrtz * zeta2) / 9. + (2. * Dilogsqrtz * zeta2) / 9. - (log2z * zeta2) / 18. - (M_PI2 * zeta2) / 27. +
                    (log_mub_Mb * (12. * Dilogminussqrtz + 12. * Dilogsqrtz + 118.990000008 * sqrtz + 6. * Dilogminussqrtz * sqrtz - 6. * Dilogsqrtz * sqrtz + 3. * M_PI2 * sqrtz - 9.624000006 * nl * sqrtz + 13.037999994 * sqrtz3 - 6. * Dilogminussqrtz * sqrtz3 +
                                   6. * Dilogsqrtz * sqrtz3 - 3. * M_PI2 * sqrtz3 - 1.206 * nl * sqrtz3 - 38.372000004 * z + 3.96 * nl * z - 3. * logz * (-1. + z) * (log1minussqrtz * (2. - sqrtz + 2. * z) + log1plussqrtz * (2. + sqrtz + 2. * z)) - 12. * Dilogminussqrtz * zeta2 -
                                   12. * Dilogsqrtz * zeta2 + 3. * log2z * zeta2 + 2. * M_PI2 * zeta2)) /
                        18. +
                    (logz * (log1minussqrtz * oneminussqrtz2 * (1. + sqrtz + z) + 1. * (sqrtz * (-76.2645000015 + 4.9559999985 * nl - 13.544000001 * sqrtz + 0.1544999985 * z) +
                                                                                        log1plussqrtz * (1. + sqrtz + sqrtz3 + zeta2)))) /
                        9. +
                    (nl * (-246643. + 288. * log2_4 + 36. * (-967. - 88. * log2 + 16. * log2_2) * M_PI2 + 732. * M_PI4 + 6912. * polylog4_12 - 78084. * zeta3)) / 23328. + (58. * zeta3) / 27. - (1439. * M_PI2 * zeta3) / 432. +
                    (nl2 * (2353. + 936. * M_PI2 + 3024. * zeta3)) / 23328. + (log_mub_Mb * (177305. + 10656. * (3. + log2) * zeta2 + 4. * nl2 * (89. + 72. * zeta2) - 4824. * zeta3 - 2. * nl * (72. * (49. + 4. * log2) * zeta2 + 7. * (1447. + 144. * zeta3)))) / 2592. + (1975. * zeta5) / 216.) /
                   M_PI3;
    // a_s(mu_b) to a_s(mu_1)
    double beta0 = 23. / 3.;
    double beta1 = 51. / 8. - 19. / 24. * 5.;
    double log_mu1_mub = 2. * log(mu_1 / mu_b);
    double as_running1 = beta0 * log_mu1_mub;
    double as_running2 = -beta0 * beta0 * log_mu1_mub * log_mu1_mub + beta1 * log_mu1_mub;
 
    for (quarks qq = cc; qq <= uu; qq = quarks(qq + 1))
    {
        for (int i = 1; i <= 8; i++)
        {
            if (i == 7)
                i++;
            for (int j = i; j <= 8; j++)
            {
                if (j == 7)
                    j++;
                if (j < i)
                    continue;
                cache_p[index_p(qq, i, j, 3)] = 8. * PoletoMS_as1 * p(qq, i, j, 2) + (32. * PoletoMS_as2 + 16. * PoletoMS_as1 * PoletoMS_as1 + 8. * PoletoMS_as1 * as_running1) * p(qq, i, j, 1) + (128. * (PoletoMS_as3 + PoletoMS_as1 * PoletoMS_as2) + 8. * PoletoMS_as1 * as_running2 + 32. * PoletoMS_as2 * as_running1) * p(qq, i, j, 0);
                cache_ps[index_p(qq, i, j, 3)] = 8. * PoletoMS_as1 * p_s(qq, i, j, 2) + (32. * PoletoMS_as2 + 16. * PoletoMS_as1 * PoletoMS_as1 + 8. * PoletoMS_as1 * as_running1) * p_s(qq, i, j, 1) + (128. * (PoletoMS_as3 + PoletoMS_as1 * PoletoMS_as2) + 8. * PoletoMS_as1 * as_running2 + 32. * PoletoMS_as2 * as_running1) * p_s(qq, i, j, 0);
                if (j >= 3 and j <= 6)
                {
                    cache_p[index_p(qq, i, j, 2)] = 0.;
                    cache_ps[index_p(qq, i, j, 2)] = 0.;
                }
                cache_p[index_p(qq, i, j, 2)] += 8. * PoletoMS_as1 * p(qq, i, j, 1) + (32. * PoletoMS_as2 + 16. * PoletoMS_as1 * PoletoMS_as1 +
                                                                                       8. * PoletoMS_as1 * as_running1) *
                                                                                          p(qq, i, j, 0);
                cache_ps[index_p(qq, i, j, 2)] += 8. * PoletoMS_as1 * p_s(qq, i, j, 1) + (32. * PoletoMS_as2 + 16. * PoletoMS_as1 * PoletoMS_as1 +
                                                                                          8. * PoletoMS_as1 * as_running1) *
                                                                                             p_s(qq, i, j, 0);
                cache_p_LO[index_p(qq, i, j, 2)] = cache_p[index_p(qq, i, j, 2)];
                cache_ps_LO[index_p(qq, i, j, 2)] = cache_ps[index_p(qq, i, j, 2)];
 
                cache_p[index_p(qq, i, j, 1)] += 8. * PoletoMS_as1 * p(qq, i, j, 0);
                cache_ps[index_p(qq, i, j, 1)] += 8. * PoletoMS_as1 * p_s(qq, i, j, 0);
                cache_p_LO[index_p(qq, i, j, 1)] = cache_p[index_p(qq, i, j, 1)];
                cache_ps_LO[index_p(qq, i, j, 1)] = cache_ps[index_p(qq, i, j, 1)];
                cache_p_LO[index_p(qq, i, j, 0)] = cache_p[index_p(qq, i, j, 0)];
                cache_ps_LO[index_p(qq, i, j, 0)] = cache_ps[index_p(qq, i, j, 0)];
            }
        }
    }
    return;
}

◆ poletoPS_pp_s()

void AmpDB2::poletoPS_pp_s ( )

private

Definition at line 2372 of file AmpDB2.cpp.

{
    // analog to arxiv:2106.05979 eq. (33) for PS mass
    double log_mu1_Mb = 2. * log(mu_1 / Mb);
    double log_mub_Mb = 2. * log(mu_b / Mb);
    // hep-ph/9804241v2 eq. (21)
    PoletoPS_as1 = 4. / 3. * mu_f / Mb;
    PoletoPS_as2 = 4. / 3. * mu_f / (Mb * 4.) * (a1 - b0 * (2. * log(mu_f / Mb) - 2.));
    bool flag_LOz = true;
    double beta0 = 23. / 3.;
    double as_running1 = beta0 * log_mu1_Mb; // a_s(Mb) to a_s(mu_1)
    PoletoMS_as1 = 4. / 3. + log_mub_Mb;
    for (quarks qq = cc; qq <= uu; qq = quarks(qq + 1))
    {
        for (int i = 1; i <= 6; i++)
        {
            // not all terms used for n=2
            for (int j = i; j <= 8; j++)
            {
                if (j == 3)
                    j = 8;
                if (i >= 3)
                    j = 8;
                cache_p[index_p(qq, i, j, 2)] += 8. * PoletoPS_as1 * p(qq, i, j, 1, flag_LOz) + (32. * PoletoPS_as2 + 16. * PoletoPS_as1 * PoletoPS_as1 +
                                                                                                 8. * PoletoPS_as1 * as_running1 -
                                                                                                 8. * PoletoPS_as1 * 4. * (-PoletoPS_as1 + PoletoMS_as1) // internal Mb_PS to Mb(mu_b)
                                                                                                 ) * p(qq, i, j, 0, flag_LOz);
                cache_ps[index_p(qq, i, j, 2)] += 8. * PoletoPS_as1 * p_s(qq, i, j, 1, flag_LOz) + (32. * PoletoPS_as2 + 16. * PoletoPS_as1 * PoletoPS_as1 +
                                                                                                    8. * PoletoPS_as1 * as_running1 -
                                                                                                    8. * PoletoPS_as1 * 4. * (-PoletoPS_as1 + PoletoMS_as1)) *
                                                                                                       p_s(qq, i, j, 0, flag_LOz);
                cache_p_LO[index_p(qq, i, j, 2)] = cache_p[index_p(qq, i, j, 2)];
                cache_ps_LO[index_p(qq, i, j, 2)] = cache_ps[index_p(qq, i, j, 2)];
            }
            for (int j = i; j <= 8; j++)
            {
                cache_p_LO[index_p(qq, i, j, 1)] += 8. * PoletoPS_as1 * p(qq, i, j, 0, flag_LOz);
                cache_ps_LO[index_p(qq, i, j, 1)] += 8. * PoletoPS_as1 * p_s(qq, i, j, 0, flag_LOz);
                if (j == 1 or j == 2 or j == 8)
                {
                    cache_p[index_p(qq, i, j, 1)] += 8. * PoletoPS_as1 * p(qq, i, j, 0);
                    cache_ps[index_p(qq, i, j, 1)] += 8. * PoletoPS_as1 * p_s(qq, i, j, 0);
                }
                else if (3 <= j and j <= 6)
                {
                    cache_p[index_p(qq, i, j, 1)] += 8. * PoletoPS_as1 * p(qq, i, j, 0, flag_LOz);
                    cache_ps[index_p(qq, i, j, 1)] += 8. * PoletoPS_as1 * p_s(qq, i, j, 0, flag_LOz);
                }
            }
        }
    }
    return;
}

◆ poletoPS_pp_s_partialN3LO()

void AmpDB2::poletoPS_pp_s_partialN3LO ( )

private

Definition at line 2502 of file AmpDB2.cpp.

{
    // analog to arxiv:2106.05979 eq. (33) for PS mass
    double log_mu1_Mb = 2. * log(mu_1 / Mb);
    double log_mub_Mb = 2. * log(mu_b / Mb);
    // hep-ph/9804241v2 eq. (21)
    double log_muf_Mb = 2. * log(mu_f / Mb);
    PoletoPS_as1 = 4. / 3. * mu_f / Mb;
    double PoletoPS_as1_2 = PoletoPS_as1 * PoletoPS_as1;
    PoletoPS_as2 = 4. / 3. * mu_f / (Mb * 4.) * (a1 - b0 * (log_muf_Mb - 2.));
    PoletoPS_as3 = 4. / 3. * mu_f / (Mb * 16.) * (a2 - (2. * a1 * b0 + b1) * (log_muf_Mb - 2.) + b0_2 * (log_muf_Mb * log_muf_Mb - 4. * log_muf_Mb + 8.));
    PoletoMS_as1 = 4. / 3. + log_mub_Mb;
    PoletoMS_as2 = (307. / 32. + M_PI2 / 3. + M_PI2 / 9. * log2 - 1. / 6. * zeta3 - (71. / 144. + M_PI2 / 18.) * nl + log_mub_Mb * (493. / 72. - nl * 13. / 36. + log_mub_Mb * (43. / 24. - nl / 12.))) + 4. / 3. * M_PI2 / 8. * sqrtz - z;
    double beta0 = 23. / 3.;
    double beta1 = 51. / 8. - 19. / 24. * 5.;
    double as_running1 = beta0 * log_mu1_Mb;
    double as_running1_mub = beta0 * 2. * log(mu_1 / mu_b);
    double as_running2 = -beta0 * beta0 * log_mu1_Mb * log_mu1_Mb + beta1 * log_mu1_Mb;
    for (quarks qq = cc; qq <= uu; qq = quarks(qq + 1))
    {
        for (int i = 1; i <= 8; i++)
        {
            if (i == 7)
                i++;
            for (int j = i; j <= 8; j++)
            {
                if (j == 7)
                    j++;
                if (j < i)
                    continue;
                cache_p[index_p(qq, i, j, 3)] = 8. * PoletoPS_as1 * p(qq, i, j, 2) + (32. * PoletoPS_as2 + 16. * PoletoPS_as1_2 + 8. * PoletoPS_as1 * as_running1 - 8. * PoletoPS_as1 * 4. * (-PoletoPS_as1 + PoletoMS_as1)) * p(qq, i, j, 1) + (128. * (PoletoPS_as3 + PoletoPS_as1 * PoletoPS_as2) - (-PoletoPS_as1 + PoletoMS_as1) * (128. * (PoletoPS_as1_2 + PoletoPS_as2) + 32. * PoletoPS_as1 * (as_running1 + as_running1_mub)) - (-PoletoPS_as2 + PoletoMS_as2) * 128. * PoletoPS_as1 + (32. * PoletoPS_as1_2 + 64. * PoletoPS_as2) * as_running1 + 8. * PoletoPS_as1 * as_running2) * p(qq, i, j, 0);
                cache_ps[index_p(qq, i, j, 3)] = 8. * PoletoPS_as1 * p_s(qq, i, j, 2) + (32. * PoletoPS_as2 + 16. * PoletoPS_as1_2 + 8. * PoletoPS_as1 * as_running1 - 8. * PoletoPS_as1 * 4. * (-PoletoPS_as1 + PoletoMS_as1)) * p_s(qq, i, j, 1) + (128. * (PoletoPS_as3 + PoletoPS_as1 * PoletoPS_as2) - (-PoletoPS_as1 + PoletoMS_as1) * (128. * (PoletoPS_as1_2 + PoletoPS_as2) + 32. * PoletoPS_as1 * (as_running1 + as_running1_mub)) - (-PoletoPS_as2 + PoletoMS_as2) * 128. * PoletoPS_as1 + (32. * PoletoPS_as1_2 + 64. * PoletoPS_as2) * as_running1 + 8. * PoletoPS_as1 * as_running2) * p_s(qq, i, j, 0);
                if (j >= 3 and j <= 6)
                {
                    cache_p[index_p(qq, i, j, 2)] = 0.;
                    cache_ps[index_p(qq, i, j, 2)] = 0.;
                }
                cache_p[index_p(qq, i, j, 2)] += 8. * PoletoPS_as1 * p(qq, i, j, 1) + (32. * PoletoPS_as2 + 16. * PoletoPS_as1_2 +
                                                                                       8. * PoletoPS_as1 * as_running1 -
                                                                                       8. * PoletoPS_as1 * 4. * (-PoletoPS_as1 + PoletoMS_as1)) *
                                                                                          p(qq, i, j, 0);
                cache_ps[index_p(qq, i, j, 2)] += 8. * PoletoPS_as1 * p_s(qq, i, j, 1) + (32. * PoletoPS_as2 + 16. * PoletoPS_as1_2 +
                                                                                          8. * PoletoPS_as1 * as_running1 -
                                                                                          8. * PoletoPS_as1 * 4. * (-PoletoPS_as1 + PoletoMS_as1)) *
                                                                                             p_s(qq, i, j, 0);
                cache_p_LO[index_p(qq, i, j, 2)] = cache_p[index_p(qq, i, j, 2)];
                cache_ps_LO[index_p(qq, i, j, 2)] = cache_ps[index_p(qq, i, j, 2)];
 
                cache_p[index_p(qq, i, j, 1)] += 8. * PoletoPS_as1 * p(qq, i, j, 0);
                cache_ps[index_p(qq, i, j, 1)] += 8. * PoletoPS_as1 * p_s(qq, i, j, 0);
                cache_p_LO[index_p(qq, i, j, 1)] = cache_p[index_p(qq, i, j, 1)];
                cache_ps_LO[index_p(qq, i, j, 1)] = cache_ps[index_p(qq, i, j, 1)];
                cache_p_LO[index_p(qq, i, j, 0)] = cache_p[index_p(qq, i, j, 0)];
                cache_ps_LO[index_p(qq, i, j, 0)] = cache_ps[index_p(qq, i, j, 0)];
            }
        }
    }
    return;
}

◆ RBd()

gslpp::complex AmpDB2::RBd ( orders order )

protected

A method to compute the ratio of the absolute value of the $B_d$ mixing amplitude over the Standard Model value.

Parameters

[in] order the QCD order of the computation

Returns: $\vert (M_{21}^{bd})_\mathrm{full}/(M_{21}^{bd})_\mathrm{SM}\vert$

Definition at line 110 of file AmpDB2.cpp.

{
    mySM.getFlavour().getHDF2().getCoeffBd().getOrder();
    if (mySM.getFlavour().getHDF2().getCoeffBd().getOrder() < getHighest(order))
        throw std::runtime_error("DmBd::computeThValue(): requires cofficient of order not computed");
 
    gslpp::vector<gslpp::complex> **allcoeff_SM = mySM.getFlavour().ComputeCoeffBd(
        mySM.getBBd().getMu(),
        mySM.getBBd().getScheme(), true);
 
    C_1_SM = ((*(allcoeff_SM[LO]))(0) + (*(allcoeff_SM[NLO]))(0));
 
    gslpp::vector<gslpp::complex> **allcoeff = mySM.getFlavour().ComputeCoeffBd(
        mySM.getBBd().getMu(),
        mySM.getBBd().getScheme());
 
    gslpp::vector<double> me(mySM.getBBd().getBpars());
    double MBd = mySM.getMesons(QCD::B_D).getMass();
    double Mb = mySM.Mrun(mySM.getBBd().getMu(),
                          mySM.getQuarks(QCD::BOTTOM).getMass_scale(),
                          mySM.getQuarks(QCD::BOTTOM).getMass(), QCD::BOTTOM, FULLNNLO);
    double Md = mySM.Mrun(mySM.getBBd().getMu(),
                          mySM.getQuarks(QCD::DOWN).getMass_scale(),
                          mySM.getQuarks(QCD::DOWN).getMass(), QCD::DOWN, FULLNNLO);
    double KBd = MBd / (Mb + Md) * MBd / (Mb + Md);
    double Fbd = mySM.getMesons(QCD::B_D).getDecayconst();
    me(0) *= 1. / 3. * MBd * Fbd * Fbd;
    me(1) *= -5. / 24. * KBd * MBd * Fbd * Fbd;
    me(2) *= 1. / 24. * KBd * MBd * Fbd * Fbd;
    me(3) *= 1. / 4. * KBd * MBd * Fbd * Fbd;
    me(4) *= 1. / 12. * KBd * MBd * Fbd * Fbd;
 
    /*std::cout << "low scale :" << std::endl << std::endl;
    std::cout << "C1_SM :" << C_1_SM << std::endl << std::endl;
 
    std::cout << "C1 :" << ((*(allcoeff[LO]))(0) + (*(allcoeff[NLO]))(0)) << std::endl;
    std::cout << "C2 :" << ((*(allcoeff[LO]))(1) + (*(allcoeff[NLO]))(1)) << std::endl;
    std::cout << "C3 :" << ((*(allcoeff[LO]))(2) + (*(allcoeff[NLO]))(2)) << std::endl;
    std::cout << "C4 :" << ((*(allcoeff[LO]))(3) + (*(allcoeff[NLO]))(3)) << std::endl;
    std::cout << "C5 :" << ((*(allcoeff[LO]))(4) + (*(allcoeff[NLO]))(4)) << std::endl << std::endl;*/
 
    switch (order)
    {
    case FULLNLO:
        return (*(allcoeff[LO]) + *(allcoeff[NLO])) * me /
               (C_1_SM * me(0));
    case LO:
        return ((*(allcoeff[LO])) * me / ((*(allcoeff_SM[LO]))(0) * me(0)));
    default:
        throw std::runtime_error("RBd::RBd(): order not implemented");
    }
}

◆ RBs()

gslpp::complex AmpDB2::RBs ( orders order )

protected

A method to compute the ratio of the absolute value of the $B_s$ mixing amplitude over the Standard Model value.

Parameters

[in] order the QCD order of the computation

Returns: $\vert (M_{21}^{bs})_\mathrm{full}/(M_{21}^{bs})_\mathrm{SM}\vert$

Definition at line 57 of file AmpDB2.cpp.

{
    mySM.getFlavour().getHDF2().getCoeffBs().getOrder();
    if (mySM.getFlavour().getHDF2().getCoeffBs().getOrder() < getHighest(order))
        throw std::runtime_error("DmBd::computeThValue(): requires cofficient of order not computed");
 
    gslpp::vector<gslpp::complex> **allcoeff_SM = mySM.getFlavour().ComputeCoeffBs(
        mySM.getBBs().getMu(),
        mySM.getBBs().getScheme(), true);
 
    C_1_SM = ((*(allcoeff_SM[LO]))(0) + (*(allcoeff_SM[NLO]))(0));
 
    gslpp::vector<gslpp::complex> **allcoeff = mySM.getFlavour().ComputeCoeffBs(
        mySM.getBBs().getMu(),
        mySM.getBBs().getScheme());
 
    gslpp::vector<double> me(mySM.getBBs().getBpars());
    double MBs = mySM.getMesons(QCD::B_S).getMass();
    double Mb = mySM.Mrun(mySM.getBBs().getMu(),
                          mySM.getQuarks(QCD::BOTTOM).getMass_scale(),
                          mySM.getQuarks(QCD::BOTTOM).getMass(), QCD::BOTTOM, FULLNNLO);
    double Ms = mySM.Mrun(mySM.getBBs().getMu(),
                          mySM.getQuarks(QCD::STRANGE).getMass_scale(),
                          mySM.getQuarks(QCD::STRANGE).getMass(), QCD::STRANGE,FULLNNLO);
    double KBs = MBs / (Mb + Ms) * MBs / (Mb + Ms);
    double Fbs = mySM.getMesons(QCD::B_S).getDecayconst();
    me(0) *= 1. / 3. * MBs * Fbs * Fbs;
    me(1) *= -5. / 24. * KBs * MBs * Fbs * Fbs;
    me(2) *= 1. / 24. * KBs * MBs * Fbs * Fbs;
    me(3) *= 1. / 4. * KBs * MBs * Fbs * Fbs;
    me(4) *= 1. / 12. * KBs * MBs * Fbs * Fbs;
 
    /*std::cout << "low scale :" << std::endl << std::endl;
    std::cout << "C1_SM :" << C_1_SM << std::endl << std::endl;
 
    std::cout << "C1 :" << ((*(allcoeff[LO]))(0) + (*(allcoeff[NLO]))(0)) << std::endl;
    std::cout << "C2 :" << ((*(allcoeff[LO]))(1) + (*(allcoeff[NLO]))(1)) << std::endl;
    std::cout << "C3 :" << ((*(allcoeff[LO]))(2) + (*(allcoeff[NLO]))(2)) << std::endl;
    std::cout << "C4 :" << ((*(allcoeff[LO]))(3) + (*(allcoeff[NLO]))(3)) << std::endl;
    std::cout << "C5 :" << ((*(allcoeff[LO]))(4) + (*(allcoeff[NLO]))(4)) << std::endl << std::endl;*/
 
    switch (order)
    {
    case FULLNLO:
        return (*(allcoeff[LO]) + *(allcoeff[NLO])) * me /
               (C_1_SM * me(0));
    case LO:
        return ((*(allcoeff[LO])) * me /(*(allcoeff_SM[LO]))(0) * me(0));
    default:
        throw std::runtime_error("RBs::RBs(): order not implemented");
    }
}

◆ transformation()

gslpp::vector< gslpp::complex > AmpDB2::transformation	(	gslpp::vector< gslpp::complex >	result,
		orders	order
	)

private

Member Data Documentation

◆ a1

double AmpDB2::a1 = 31./3. - 10./9. * nl

private

Definition at line 282 of file AmpDB2.h.

◆ a2

double AmpDB2::a2

private

Initial value:

= (4343./162. + 6.*M_PI2 - M_PI4/4. + 22./3. * zeta3) * 3. * 3. - (1798./81. + 56./3. * zeta3) * 3. * 0.5 * nl

- (55./3. - 16. * zeta3) * 4./3. * 0.5 * nl + 400./81. * 0.25 * nl * nl

Definition at line 283 of file AmpDB2.h.

◆ as_4pi

double AmpDB2::as_4pi

private

Definition at line 259 of file AmpDB2.h.

◆ as_4pi_mu1

double AmpDB2::as_4pi_mu1

private

Definition at line 257 of file AmpDB2.h.

◆ as_4pi_mu2

double AmpDB2::as_4pi_mu2

private

Definition at line 258 of file AmpDB2.h.

◆ b0

double AmpDB2::b0 = 11. - 2. * nl/3.

private

Definition at line 285 of file AmpDB2.h.

◆ b0_2

double AmpDB2::b0_2 = b0 * b0

private

Definition at line 286 of file AmpDB2.h.

◆ b1

double AmpDB2::b1 = 102. - 38. * nl/3.

private

Definition at line 287 of file AmpDB2.h.

◆ BMeson

int AmpDB2::BMeson

private

B meson index (0 for Bd and 1 for Bs)

Definition at line 171 of file AmpDB2.h.

◆ C1_LO_1overm

gslpp::complex AmpDB2::C1_LO_1overm

private

Definition at line 467 of file AmpDB2.h.

◆ C2_LO_1overm

gslpp::complex AmpDB2::C2_LO_1overm

private

Definition at line 468 of file AmpDB2.h.

◆ C_1_SM

gslpp::complex AmpDB2::C_1_SM

private

Wilson coeffients H_{eff}^{DF2} $C_1$

Definition at line 173 of file AmpDB2.h.

◆ C_Buras_LO

gslpp::complex AmpDB2::C_Buras_LO[8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

private

Definition at line 302 of file AmpDB2.h.

◆ C_Buras_NLO

gslpp::complex AmpDB2::C_Buras_NLO[8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

private

Definition at line 303 of file AmpDB2.h.

◆ C_Buras_NNLO

gslpp::complex AmpDB2::C_Buras_NNLO[8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

private

Definition at line 304 of file AmpDB2.h.

◆ C_Misiak_LO

gslpp::complex AmpDB2::C_Misiak_LO[8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

private

Definition at line 299 of file AmpDB2.h.

◆ C_Misiak_NLO

gslpp::complex AmpDB2::C_Misiak_NLO[8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

private

Definition at line 300 of file AmpDB2.h.

◆ C_Misiak_NNLO

gslpp::complex AmpDB2::C_Misiak_NNLO[8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

private

Definition at line 301 of file AmpDB2.h.

◆ cache_deltas_1overm

gslpp::complex AmpDB2::cache_deltas_1overm[6] = {0.}

private

Definition at line 452 of file AmpDB2.h.

◆ cache_deltas_1overm_tradBasis

gslpp::complex AmpDB2::cache_deltas_1overm_tradBasis[6] = {0.}

private

Definition at line 364 of file AmpDB2.h.

◆ cache_p

double AmpDB2::cache_p[768] = { 0. }

private

Definition at line 400 of file AmpDB2.h.

◆ cache_p_LO

double AmpDB2::cache_p_LO[576] = {0.}

private

Definition at line 404 of file AmpDB2.h.

◆ cache_ps

double AmpDB2::cache_ps[768] = { 0. }

private

Definition at line 401 of file AmpDB2.h.

◆ cache_ps_LO

double AmpDB2::cache_ps_LO[576] = {0.}

private

Definition at line 405 of file AmpDB2.h.

◆ cacheC

gslpp::complex AmpDB2::cacheC[8] = { 0., 0., 0., 0., 0., 0., NAN, 0.}

private

Definition at line 298 of file AmpDB2.h.

◆ cacheD

gslpp::complex AmpDB2::cacheD[6] = {0.}

private

Definition at line 332 of file AmpDB2.h.

◆ cacheF0

double AmpDB2::cacheF0[24] = {0.}

private

Definition at line 329 of file AmpDB2.h.

◆ cacheF1

double AmpDB2::cacheF1[24] = {0.}

private

Definition at line 330 of file AmpDB2.h.

◆ cacheg

gslpp::complex AmpDB2::cacheg[12] = {0.}

private

Definition at line 460 of file AmpDB2.h.

◆ cachegtilde

gslpp::complex AmpDB2::cachegtilde[12] = {0.}

private

Definition at line 461 of file AmpDB2.h.

◆ cacheP

double AmpDB2::cacheP[84] = {0.}

private

Definition at line 331 of file AmpDB2.h.

◆ Cl2PI3

const double AmpDB2::Cl2PI3 = 1.014941606

private

Definition at line 193 of file AmpDB2.h.

◆ coeffsMStoRI

gslpp::matrix<double> AmpDB2::coeffsMStoRI

private

Definition at line 212 of file AmpDB2.h.

◆ Dilogsigma

double AmpDB2::Dilogsigma

private

Definition at line 255 of file AmpDB2.h.

◆ Dilogsigma2

double AmpDB2::Dilogsigma2

private

Definition at line 256 of file AmpDB2.h.

◆ Dilogz

double AmpDB2::Dilogz

private

Definition at line 254 of file AmpDB2.h.

◆ flag_fixmub

bool AmpDB2::flag_fixmub

private

Definition at line 213 of file AmpDB2.h.

◆ flag_LOz

bool AmpDB2::flag_LOz = true

private

Definition at line 403 of file AmpDB2.h.

◆ flag_resumz

bool AmpDB2::flag_resumz

private

Definition at line 207 of file AmpDB2.h.

◆ flag_RI

bool AmpDB2::flag_RI

private

Definition at line 214 of file AmpDB2.h.

◆ Gf2

double AmpDB2::Gf2

private

Definition at line 235 of file AmpDB2.h.

◆ K_1

gslpp::complex AmpDB2::K_1

private

Definition at line 470 of file AmpDB2.h.

◆ K_2

gslpp::complex AmpDB2::K_2

private

Definition at line 471 of file AmpDB2.h.

◆ lambda_c_d

gslpp::complex AmpDB2::lambda_c_d

private

Definition at line 374 of file AmpDB2.h.

◆ lambda_c_s

gslpp::complex AmpDB2::lambda_c_s

private

Definition at line 376 of file AmpDB2.h.

◆ lambda_u_d

gslpp::complex AmpDB2::lambda_u_d

private

Definition at line 375 of file AmpDB2.h.

◆ lambda_u_s

gslpp::complex AmpDB2::lambda_u_s

private

Definition at line 377 of file AmpDB2.h.

◆ log12sqrt52

const double AmpDB2::log12sqrt52 = log(0.5 + sqrt5 / 2.)

private

Definition at line 191 of file AmpDB2.h.

◆ log1minus4z

double AmpDB2::log1minus4z

private

Definition at line 244 of file AmpDB2.h.

◆ log1minusz

double AmpDB2::log1minusz

private

Definition at line 243 of file AmpDB2.h.

◆ log2

const double AmpDB2::log2 = log(2)

private

Definition at line 185 of file AmpDB2.h.

◆ log2_2

const double AmpDB2::log2_2 = log2 * log2

private

Definition at line 186 of file AmpDB2.h.

◆ log2_4

const double AmpDB2::log2_4 = log2_2 * log2_2

private

Definition at line 187 of file AmpDB2.h.

◆ log2sigma

double AmpDB2::log2sigma

private

Definition at line 249 of file AmpDB2.h.

◆ log2z

double AmpDB2::log2z

private

Definition at line 242 of file AmpDB2.h.

◆ log3

const double AmpDB2::log3 = log(3)

private

Definition at line 188 of file AmpDB2.h.

◆ logsigma

double AmpDB2::logsigma

private

Definition at line 248 of file AmpDB2.h.

◆ logx_1

double AmpDB2::logx_1

private

Definition at line 252 of file AmpDB2.h.

◆ logx_2

double AmpDB2::logx_2

private

Definition at line 253 of file AmpDB2.h.

◆ logz

double AmpDB2::logz

private

Definition at line 241 of file AmpDB2.h.

◆ M_PI2

const double AmpDB2::M_PI2 = M_PI * M_PI

private

Definition at line 178 of file AmpDB2.h.

◆ M_PI3

const double AmpDB2::M_PI3 = M_PI2 * M_PI

private

Definition at line 179 of file AmpDB2.h.

◆ M_PI4

const double AmpDB2::M_PI4 = M_PI2 * M_PI2

private

Definition at line 180 of file AmpDB2.h.

◆ Mb

double AmpDB2::Mb

private

Definition at line 270 of file AmpDB2.h.

◆ MB

double AmpDB2::MB

private

Definition at line 271 of file AmpDB2.h.

◆ Mb2_prefactor

double AmpDB2::Mb2_prefactor

private

Definition at line 273 of file AmpDB2.h.

◆ Mb2_prefactor_1overm

double AmpDB2::Mb2_prefactor_1overm

private

Definition at line 274 of file AmpDB2.h.

◆ Mb_Mb

double AmpDB2::Mb_Mb

private

Definition at line 275 of file AmpDB2.h.

◆ Mb_pole

double AmpDB2::Mb_pole

private

Definition at line 276 of file AmpDB2.h.

◆ Mb_PS

double AmpDB2::Mb_PS

private

Definition at line 277 of file AmpDB2.h.

◆ MB_s

double AmpDB2::MB_s

private

Definition at line 272 of file AmpDB2.h.

◆ Mc

double AmpDB2::Mc

private

Definition at line 269 of file AmpDB2.h.

◆ Md

double AmpDB2::Md

private

Definition at line 267 of file AmpDB2.h.

◆ me

gslpp::vector<double> AmpDB2::me = gslpp::vector<double>(5, 0.)

private

Definition at line 201 of file AmpDB2.h.

◆ me_R

gslpp::vector<double> AmpDB2::me_R = gslpp::vector<double>(5, 0.)

private

Definition at line 203 of file AmpDB2.h.

◆ me_Rtilde

gslpp::vector<double> AmpDB2::me_Rtilde = gslpp::vector<double>(3, 0.)

private

Definition at line 204 of file AmpDB2.h.

◆ meMStoRI

gslpp::matrix<double> AmpDB2::meMStoRI

private

Definition at line 210 of file AmpDB2.h.

◆ meoverme0

gslpp::vector<double> AmpDB2::meoverme0 = gslpp::vector<double>(3, 0.)

private

Definition at line 202 of file AmpDB2.h.

◆ Ms

double AmpDB2::Ms

private

Definition at line 268 of file AmpDB2.h.

◆ mu_1

double AmpDB2::mu_1

private

Definition at line 196 of file AmpDB2.h.

◆ mu_1_overm

double AmpDB2::mu_1_overm

private

Definition at line 197 of file AmpDB2.h.

◆ mu_2

double AmpDB2::mu_2

private

Definition at line 198 of file AmpDB2.h.

◆ mu_b

double AmpDB2::mu_b

private

Definition at line 199 of file AmpDB2.h.

◆ mu_f

double AmpDB2::mu_f = 2.

private

Definition at line 280 of file AmpDB2.h.

◆ mySM

const StandardModel& AmpDB2::mySM

private

Model type

Definition at line 169 of file AmpDB2.h.

◆ nl

double AmpDB2::nl = 4.

private

Definition at line 281 of file AmpDB2.h.

◆ oneminusz2

double AmpDB2::oneminusz2

private

Definition at line 245 of file AmpDB2.h.

◆ oneminusz_1overm2

double AmpDB2::oneminusz_1overm2

private

Definition at line 264 of file AmpDB2.h.

◆ PoletoMS_as1

double AmpDB2::PoletoMS_as1

private

Definition at line 417 of file AmpDB2.h.

◆ PoletoMS_as2

double AmpDB2::PoletoMS_as2

private

Definition at line 418 of file AmpDB2.h.

◆ PoletoMS_as2_z0

double AmpDB2::PoletoMS_as2_z0

private

Definition at line 420 of file AmpDB2.h.

◆ PoletoMS_as2_z1

double AmpDB2::PoletoMS_as2_z1

private

Definition at line 421 of file AmpDB2.h.

◆ PoletoMS_as3

double AmpDB2::PoletoMS_as3

private

Definition at line 419 of file AmpDB2.h.

◆ PoletoPS_as1

double AmpDB2::PoletoPS_as1

private

Definition at line 427 of file AmpDB2.h.

◆ PoletoPS_as2

double AmpDB2::PoletoPS_as2

private

Definition at line 428 of file AmpDB2.h.

◆ PoletoPS_as3

double AmpDB2::PoletoPS_as3

private

Definition at line 429 of file AmpDB2.h.

◆ polylog4_12

const double AmpDB2::polylog4_12 = 0.517479062

private

Definition at line 194 of file AmpDB2.h.

◆ sigma

double AmpDB2::sigma

private

Definition at line 247 of file AmpDB2.h.

◆ sqrt1minus4z

double AmpDB2::sqrt1minus4z

private

Definition at line 246 of file AmpDB2.h.

◆ sqrt1minus4z_1overm

double AmpDB2::sqrt1minus4z_1overm

private

Definition at line 265 of file AmpDB2.h.

◆ sqrt3

const double AmpDB2::sqrt3 = sqrt(3)

private

Definition at line 189 of file AmpDB2.h.

◆ sqrt5

const double AmpDB2::sqrt5 = sqrt(5)

private

Definition at line 190 of file AmpDB2.h.

◆ sqrtz

double AmpDB2::sqrtz

private

Definition at line 240 of file AmpDB2.h.

◆ t_2

const double AmpDB2::t_2 = -0.389011713

private

Definition at line 192 of file AmpDB2.h.

◆ VcbVcd

gslpp::complex AmpDB2::VcbVcd

private

Definition at line 317 of file AmpDB2.h.

◆ VcbVcd2

gslpp::complex AmpDB2::VcbVcd2

private

Definition at line 319 of file AmpDB2.h.

◆ VcbVcs

gslpp::complex AmpDB2::VcbVcs

private

Definition at line 318 of file AmpDB2.h.

◆ VcbVcs2

gslpp::complex AmpDB2::VcbVcs2

private

Definition at line 320 of file AmpDB2.h.

◆ VtbVtd

gslpp::complex AmpDB2::VtbVtd

private

Definition at line 313 of file AmpDB2.h.

◆ VtbVtd2

gslpp::complex AmpDB2::VtbVtd2

private

Definition at line 315 of file AmpDB2.h.

◆ VtbVts

gslpp::complex AmpDB2::VtbVts

private

Definition at line 314 of file AmpDB2.h.

◆ VtbVts2

gslpp::complex AmpDB2::VtbVts2

private

Definition at line 316 of file AmpDB2.h.

◆ x_1

double AmpDB2::x_1

private

Definition at line 250 of file AmpDB2.h.

◆ x_2

double AmpDB2::x_2

private

Definition at line 251 of file AmpDB2.h.

◆ z

double AmpDB2::z

private

Definition at line 236 of file AmpDB2.h.

◆ z2

double AmpDB2::z2

private

Definition at line 237 of file AmpDB2.h.

◆ z3

double AmpDB2::z3

private

Definition at line 238 of file AmpDB2.h.

◆ z4

double AmpDB2::z4

private

Definition at line 239 of file AmpDB2.h.

◆ z_1overm

double AmpDB2::z_1overm

private

Definition at line 262 of file AmpDB2.h.

◆ z_1overm2

double AmpDB2::z_1overm2

private

Definition at line 263 of file AmpDB2.h.

◆ zeta2

const double AmpDB2::zeta2 = gslpp_special_functions::zeta(2)

private

Definition at line 181 of file AmpDB2.h.

◆ zeta3

const double AmpDB2::zeta3 = gslpp_special_functions::zeta(3)

private

Definition at line 182 of file AmpDB2.h.

◆ zeta4

const double AmpDB2::zeta4 = gslpp_special_functions::zeta(4)

private

Definition at line 183 of file AmpDB2.h.

◆ zeta5

const double AmpDB2::zeta5 = gslpp_special_functions::zeta(5)

private

Definition at line 184 of file AmpDB2.h.

The documentation for this class was generated from the following files:

Detailed Description

Public Types

Public Member Functions

Protected Member Functions

Private Member Functions

Private Attributes

Member Enumeration Documentation

◆ mass_schemes

◆ quark

◆ quarks

Constructor & Destructor Documentation

◆ AmpDB2()

Member Function Documentation

◆ c()

◆ c_H()

◆ c_H_partial()

◆ compute_deltas_1overm()

◆ compute_deltas_1overm_tradBasis()

◆ compute_g()

◆ compute_matrixelements()

◆ compute_partialNNLO()

◆ compute_pp_s()

◆ computeCKMandMasses()

◆ computeD()

◆ computeF0()

◆ computeF1()

◆ computeP()

◆ computeWilsonCoeffsBuras()

◆ computeWilsonCoeffsMisiak()

◆ D()

◆ delta_1overm()

◆ delta_1overm_tradBasis()

◆ deltas_1overm()

◆ deltas_1overm_tradBasis()

◆ F()

◆ F0()

◆ F1()

◆ g()

◆ Gamma21overM21()

◆ Gamma21overM21_tradBasis()

◆ getGamma21overM21()

◆ getGamma21overM21_tradBasis()

◆ getM21()

◆ getPB()

◆ getRB()

◆ gtilde()

◆ H()

◆ H_allpartial()

◆ H_partial()

◆ H_s()

◆ H_s_allpartial()

◆ H_s_partial()

◆ index_deltas()

◆ index_p()

◆ indexD()

◆ indexF()

◆ indexg()

◆ indexP()

◆ M21_Bd()

◆ M21_Bs()

◆ p()

◆ P()

◆ p_s()

◆ PBd()

◆ PBs()

◆ poletoMSbar_pp_s()

◆ poletoMSbar_pp_s_partialN3LO()

◆ poletoPS_pp_s()

◆ poletoPS_pp_s_partialN3LO()

◆ RBd()

◆ RBs()

◆ transformation()

Member Data Documentation

◆ a1

◆ a2

◆ as_4pi

◆ as_4pi_mu1

◆ as_4pi_mu2

◆ b0

◆ b0_2