file FlavBit/flav_loop_functions.hpp
[No description available] More…
Namespaces
| Name |
|---|
| Gambit TODO: see if we can use this one: |
| Gambit::FlavBit |
| Gambit::FlavBit::LoopFunctions |
| Gambit::FlavBit::Vertices |
| Gambit::FlavBit::Penguins |
| Gambit::FlavBit::Boxes |
| Gambit::FlavBit::FormFactors |
Detailed Description
Author: Tomas Gonzalo (t.e.gonzalo@fys.uio.no)
Date: 2017 Aug, 2018 Feb
Loop functions for flavour violating decays of charged leptons (from hep-ph/9403398) And for RK from 1706.07570
Authors (add name and date if you modify):
Source code
// GAMBIT: Global and Modular BSM Inference Tool
// *********************************************
/// \file
///
/// Loop functions for flavour violating decays of charged leptons (from hep-ph/9403398)
/// And for RK from 1706.07570
///
/// *********************************************
///
/// Authors (add name and date if you modify):
///
/// \author Tomas Gonzalo
/// (t.e.gonzalo@fys.uio.no)
/// \date 2017 Aug, 2018 Feb
///
/// *********************************************
#ifndef __flav_loop_functions_hpp
#define __flav_loop_functions_hpp
namespace Gambit
{
namespace FlavBit
{
// Loop functions for LFV diagrams
namespace LoopFunctions
{
double G1(const double x)
{
if(x == 0)
return -7./12.;
if(x == 1)
return -5./12.;
else
return (-7. + 33.*x - 57.*x*x + 31.*x*x*x + 6.*x*x*(1. - 3*x)*log(x))/(12.*pow(1.-x,4));
}
double G1(const double a, const double b, const double c)
{
if(b == c and b != 0)
return G1(a/b)/b;
else
return 0; // TODO: 2C12 + 2C22 - C1 or 2C12 + C11 - C2
}
double MFVV(const double a, const double b)
{
if(a == b)
return 1. / (3. * b);
else if(a == 0)
return 5. / (9. * b);
else
return (6.*a*a*(a-3.*b)*log(a/b) - (a-b)*(5.*a*a - 22.*a*b + 5.*b*b))/(9.*pow(a-b,4));
}
double B1(const double a, const double b, const double Q)
{
if(a == b)
return 0.5 * log(b / pow(Q,2));
else if(a == 0)
return -0.25 + 0.5*log(b / pow(Q,2));
else
return -0.5 + 0.5*log(b / pow(Q,2)) - (a*a - b*b + 2.*a*a*log(b/a)) / (4.*pow(a-b,2));
}
double B0(const double a, const double b, const double Q)
{
// TODO: behaviour when a = 0 and b = 0 undefined
if(a == 0 and b == 0)
return 0;
else if(a == b)
return -log(b / pow(Q,2));
else if(a == 0)
return 1. - log(b / pow(Q,2));
else if(b == 0)
return 1. - log(a) - log(1./pow(Q,2));
else
return 1. - log(b / pow(Q,2)) + 1./(b-a) * a * log(a/b);
}
double C0(const double a, const double b, const double c)
{
// TODO: behaviour when two paramers are 0 undefined, set it to zero
if(a == 0 and b == 0 and c == 0)
return 0;
if(a == 0 and b == 0)
return 0;
else if(a == 0 and c == 0)
return 0;
else if(b == 0 and c == 0)
return 0;
else if(c == 0)
return C0(a,c,b);
else if(a == b and b == c)
return - 1./(2*c);
else if(a == b)
return (-b + c- c*log(c/b)) / pow(b-c,2);
else if(a == c and b != 0)
return C0(a,c,b);
else if(a == c and b == 0)
return -1./c;
else if(b == c and a != 0)
return (a - c + a*log(c/a)) / pow(a-c,2);
else if(b == c and a == 0)
return -1./c;
else if(a == 0)
return (-log(b) + log(c)) / (b-c);
else if(b == 0)
return log(c/a)/(a-c);
else
return -1. / (a-b)*(a-c)*(b-c)*( b*(c-a)*log(b/a) + c*(a-b)*log(c/a));
}
double C00(const double a, const double b, const double c, const double Q)
{
// TODO: behaviour when all three parameters are zero is undefined, set it to zero
if(a == 0 and b == 0 and c == 0)
return 0;
else if(b == 0 and c == 0)
return 0.125*(3. - 2.*log(a/pow(Q,2)));
else if(a == 0 and b == 0)
return 0.125*(3. - 2.*log(c) - 2.*log(1./pow(Q,2)));
else if(c == 0)
return C00(a,c,b,Q);
else if(a == b and b == c)
return -0.25*log(c/pow(Q,2));
else if(a == b)
return - (2.*c*c*log(c/b) + (b-c)*(-b + 3.*c + 2.*(b-c)*log(b/pow(Q,2))))/(8.*pow(b-c,2));
else if(a == c and b != 0)
return C00(a,c,b,Q);
else if(a == c and b == 0)
return 0.125*(1. - 2.*log(c/pow(Q,2)));
else if(b == c and a != 0)
return (2.*(2.*a-c)*c*log(c/a)-(a-c)*(-3.*a+c+2.*(a-c)*log(a/pow(Q,2))))/(8.*pow(a-c,2));
else if(b == c and a == 0)
return 0.125*(1. - 2.*log(c) - 2.*log(1./pow(Q,2)));
else if(a == 0)
return -(2.*b*log(b) - 2.*c*log(c) + (b-c)*(-3.+2.*log(1./pow(Q,2))))/(8.*(b-c));
else if(b == 0)
return (2.*c*log(c/a) - (a-c)*(-3. + 2.*log(a/pow(Q,2))))/(8.*(a-c));
else
return 1. / (8.*(a-b)*(a-c)*(b-c)) * ( (c-a)*((a-b)*(2.*log(a/pow(Q,2))-3.)*(b-c) - 2.*b*b*log(b/a)) + 2.*c*c*(b-a)*log(c/a));
}
// Finite combination of loop functions that appears in VZw10
double B02C00C0(const double a, const double b, const double c, const double Q)
{
if(a == 0 and b == 0)
return 0.25*(1.0 - 2.0*log(c) - 2.0*log(1 / pow(Q,2)));
else
return B0(a,b,Q) - 2*C00(a,b,c,Q) + C0(a,b,c)*c;
}
double D0(const double a, const double b, const double c, const double d)
{
//TODO: behaviour when two or more parameters are zero is undefined, set it to zero
if((!a and !b) or (!b and !c) or (!b and !d) or (!c and !d))
return 0;
else if(c == 0)
return D0(a,c,b,d);
else if(d == 0)
return D0(a,d,c,b);
else if(a == b and b == c and c == d)
return 1. / (6.*d*d);
else if(a == b and b == c)
return D0(a,d,c,d);
else if(a == b and b == d)
return D0(a,c,b,d);
else if(a == c and c == d and b == 0)
return 1. / (2.*c*c);
else if(a == c and c == d and b != 0)
return (-b*b + c*c + 2.*b*c*log(b/c)) / (2.*c*pow(c-b,3));
else if(b == c and c == d and a == 0)
return 1. / (2.*d*d);
else if(b == c and c == d and a != 0)
return (a*a - d*d + 2.*a*d*log(d/a)) / (2.*d*pow(a-d,3));
else if(a == d and b == c)
return (-2*c + 2*d + (c+d)*log(c/d)) / pow(c-d,3);
else if(a == d and b == 0)
return (c - d -d*log(c/d)) / (pow(c-d,2)*d);
else if(a == d)
return 1./ ((b-d)*(d-c))-(b*log(b/d))/((b-c)*pow(b-d,2))+(c*log(c/d))/((b-c)*pow(c-d,2));
else if(a == c)
return D0(a,b,d,c);
else if(a == b)
return D0(a,d,c,b);
else if(b == c)
return D0(a,d,c,b);
else if(b == d)
return D0(a,c,b,d);
else if(c == d and b == 0)
return (a - d + d*log(d/a)) / (d*pow(a-d,2));
else if(c == d and a == 0)
return (b - d + d*log(d/b)) / (d*pow(b-d,2));
else if(c == d)
return (b*pow(a-d,2)*log(b/a) - (a-b)*( (a-d)*(b-d) + (a*b-d*d)*log(d/a) )) / ((a-b)*pow(a-d,2)*pow(b-d,2));
else if(b == 0)
return log(c/a)/((a-c)*(c-d)) + log(d/a)/((a-d)*(d-c));
else if(a == 0)
return ((d-c)*log(b) + (b-d)*log(c) + (c-b)*log(d))/((b-c)*(b-d)*(c-d));
else
return -(b*log(b/a)/((b-a)*(b-c)*(b-d)) + c*log(c/a)/((c-a)*(c-b)*(c-d)) + d*log(d/a)/((d-a)*(d-b)*(d-c)));
}
double D27(const double a, const double b, const double c, const double d)
{
//TODO: behaviour when three or more parameters are zero is undefined, set it to zero
if((!a and !b and !c) or (!a and !b and !d) or (!a and !c and !d) or (!b and !c and !d))
return 0;
if(a == b and b == c and c == d)
return -1./(12.*d);
if(a == d and c == d and b == 0)
return -1. / (8.*d);
if(a == d and c == d)
return (3.*b*b - 4.*b*d + d*d - 2.*b*b*log(b/d))/(8.*pow(b-d,3));
if(b == c and c == d)
return D27(b,a,c,d);
if(a == b and b == c)
return D27(a,d,c,b);
if(a == b and b == d)
return D27(a,c,b,d);
if(a == b and c == d and b == 0)
return -1./(4.*d);
if(a == b and c == d and d == 0)
return -1./(4.*b);
if(a == b and c == d)
return (-b*b + d*d -2.*b*d*log(d/b)) / (4.*pow(b-d,3));
if(a == c and b == d)
return D27(a,c,b,d);
if(a == d and b == c)
return D27(a,b,d,c);
if(a == b and b == 0)
return log(d/c)/(4.*(c-d));
if(a == b and c == 0)
return - (b -d +d*log(d/b))/(4.*pow(b-d,2));
if(a == b and d == 0)
return D27(a,b,d,c);
if(a == b)
return 0.25*(-c*c*log(c/b)/(pow(b-c,2)*(c-d)) + (b*(d-b)/(b-c) + d*d*log(d/b)/(c-d))/pow(b-d,2));
if(a == c)
return D27(a,c,b,d);
if(a == d)
return D27(a,d,c,b);
if(b == c and a == 0)
return (-c+d+d*log(c/d))/(4.*pow(c-d,2));
if(b == c and c == 0)
return log(d/a)/(4.*(a-d));
if(b == c and d == 0)
return (a-c+a*log(c/a))/(4.*pow(a-c,2));
if(b == c)
return (c*(a-d)*(a*(c-2.*d)+c*d)*log(c/a)+(a-c)*(c*(a-d)*(c-d)+(a-c)*d*d*log(d/a)))/(4.*pow(a-c,2)*(a-d)*pow(c-d,2));
if(b == d)
return D27(a,b,d,c);
if(c == d)
return D27(a,c,d,b);
if(a == 0)
return (b*(-c+d)*log(b)+c*(b-d)*log(c)+(-b+c)*d*log(d))/(4.*(b-c)*(b-d)*(c-d));
if(b == 0)
return ((c*log(c/a))/((a - c)*(c - d)) + (d*log(d/a))/((a - d)*(-c + d)))/4.;
if(c == 0)
return D27(a,c,b,d);
if(d == 0)
return D27(a,d,c,b);
else
return -0.25*(b*b*log(b/a)/((b-a)*(b-c)*(b-d)) + c*c*log(c/a)/((c-a)*(c-b)*(c-d)) + d*d*log(d/a)/((d-a)*(d-b)*(d-c)));
}
double IC0D0(const double a, const double b, const double c, const double d)
{
return C0(a,b,c) + d*D0(a,b,c,d);
}
}
// Loop function for RK
namespace LoopFunctions
{
double E(const double x, const double y)
{
if(x == 0 or y == 0)
return 0.0;
if(x == y)
return (x*(-4.0 + 15.0*x - 12.0*pow(x,2) + pow(x,3) + 6.0*pow(x,2)*log(x)))/ (4.*pow(-1.0 + x,3));
return x*y*(-3.0/(4.0*(1.0-x)*(1.0 - y)) + ((0.25 - 3.0/(4.0*pow(-1.0 + x,2)) - 3.0/(2.0*(-1.0 + x)))*log(x))/(x - y) + ((0.25 - 3.0/(4.0*pow(-1.0 + y,2)) - 3.0/(2.0*(-1 + y)))*log(y))/(-x + y));
}
}
// Vertices for LFV diagrams
namespace Vertices
{
// Fermion-vector vertices
std::complex<double> VpL(int i, int j, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U)
{
double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
return -1. / sqrt(2) * g2 * U(i,j);
}
double EL(int i,int j, SMInputs sminputs)
{
if(i != j) return 0;
double e = sqrt(4. * pi / sminputs.alphainv);
double g1 = e * sminputs.mZ / sminputs.mW;
double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
double cw = sminputs.mW / sminputs.mZ;
double sw = sqrt(1. - cw*cw);
return 0.5 * (-g1*sw + g2*cw);
}
double ER(int i, int j, SMInputs sminputs)
{
if(i != j) return 0;
double e = sqrt(4. * pi / sminputs.alphainv);
double g1 = e * sminputs.mZ / sminputs.mW;
double cw = sminputs.mW / sminputs.mZ;
double sw = sqrt(1. - cw*cw);
return - g1*sw;
}
std::complex<double> VL(int i, int j, SMInputs sminputs)
{
double e = sqrt(4. * pi / sminputs.alphainv);
double g1 = e * sminputs.mZ / sminputs.mW;
double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
double cw = sminputs.mW / sminputs.mZ;
double sw = sqrt(1. - cw*cw);
if(i == j)
return -0.5*(g1*sw + g2*cw);
else
return 0.;
}
std::complex<double> VR(int i, int j, SMInputs sminputs)
{
double e = sqrt(4. * pi / sminputs.alphainv);
double g1 = e * sminputs.mZ / sminputs.mW;
double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
double cw = sminputs.mW / sminputs.mZ;
double sw = sqrt(1. - cw*cw);
if(i == j)
return 0.5*(g1*sw + g2*cw);
else
return 0.;
}
std::complex<double> DL(int i, int j, SMInputs sminputs)
{
double e = sqrt(4. * pi / sminputs.alphainv);
double g1 = e * sminputs.mZ / sminputs.mW;
double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
double cw = sminputs.mW / sminputs.mZ;
double sw = sqrt(1. - cw*cw);
if(i == j)
return 1./6. * (3.*g2*cw + g1*sw);
else
return 0;
}
std::complex<double> DR(int i, int j, SMInputs sminputs)
{
double e = sqrt(4. * pi / sminputs.alphainv);
double g1 = e * sminputs.mZ / sminputs.mW;
double cw = sminputs.mW / sminputs.mZ;
double sw = sqrt(1. - cw*cw);
if(i == j)
return -1./3.*g1*sw;
else
return 0.;
}
std::complex<double> UL(int i, int j, SMInputs sminputs)
{
double e = sqrt(4. * pi / sminputs.alphainv);
double g1 = e * sminputs.mZ / sminputs.mW;
double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
double cw = sminputs.mW / sminputs.mZ;
double sw = sqrt(1. - cw*cw);
if(i == j)
return -1./6. * (3.*g2*cw - g1*sw);
else
return 0;
}
std::complex<double> UR(int i, int j, SMInputs sminputs)
{
double e = sqrt(4. * pi / sminputs.alphainv);
double g1 = e * sminputs.mZ / sminputs.mW;
double cw = sminputs.mW / sminputs.mZ;
double sw = sqrt(1. - cw*cw);
if(i == j)
return 2./3.*g1*sw;
else
return 0.;
}
std::complex<double> VuL(int i, int j, SMInputs sminputs)
{
double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
Eigen::Matrix3cd VCKM;
double lambda = sminputs.CKM.lambda, A = sminputs.CKM.A;
double rhobar = sminputs.CKM.rhobar, etabar = sminputs.CKM.etabar;
std::complex<double> I = {0,1};
std::complex<double> Vub = real(rhobar + I*etabar)*sqrt(1.-A*A*pow(lambda,4))/(sqrt(1.-pow(lambda,2))*(1.- A*A*pow(lambda,4)*(rhobar+I*etabar)));
double rho = real(Vub);
double eta = imag(Vub);
VCKM << 1. - 0.5*pow(lambda,2), lambda, A*pow(lambda,3)*(rho - I*eta),
-lambda, 1. - 0.5*pow(lambda,2), A*pow(lambda,2),
A*pow(lambda,3)*(1. - eta - I*eta), -A*pow(lambda,2), 1;
return -1./sqrt(2) * g2 * VCKM(i,j);
}
// Vector vertices
double Fw(SMInputs sminputs)
{
return sqrt(4.* pi/ sminputs.alphainv);
}
double Zww(SMInputs sminputs)
{
double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
return -g2 * sminputs.mW / sminputs.mZ;
}
// Scalar vertices
double HL(int i, int j, SMInputs sminputs)
{
double vev = 1. / sqrt(sqrt(2.)*sminputs.GF);
if(i == 0 and j == 0)
return -1. / vev * sminputs.mE;
if(i == 1 and j == 1)
return -1. / vev * sminputs.mMu;
if(i == 2 and j == 2)
return -1. / vev * sminputs.mTau;
else
return 0;
}
double HR(int i, int j, SMInputs sminputs)
{
return HL(i, j , sminputs);
}
double HdL(int i, int j, SMInputs sminputs)
{
double vev = 1. / sqrt(sqrt(2.)*sminputs.GF);
if(i == 0 and j == 0)
return -1. / vev * sminputs.mD;
else if(i == 1 and j == 1)
return -1. / vev * sminputs.mS;
else if(i == 2 and j == 2)
return -1. / vev * sminputs.mBmB;
else
return 0;
}
double HdR(int i, int j, SMInputs sminputs)
{
return HdL(i, j, sminputs);
}
double HuL(int i, int j, SMInputs sminputs)
{
double vev = 1. / sqrt(sqrt(2.)*sminputs.GF);
if(i == 0 and j == 0)
return -1. / vev * sminputs.mU;
else if(i == 1 and j == 1)
return -1. / vev * sminputs.mCmC;
else if(i == 2 and j == 2)
return -1. / vev * sminputs.mT;
else
return 0;
}
double HuR(int i, int j, SMInputs sminputs)
{
return HuL(i, j, sminputs);
}
}
// Penguin contributions
namespace Penguins
{
// Fotonic penguins
std::complex<double> A1R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
std::complex<double> a1r = {0,0};
for(int a=0; a<6; a++)
{
a1r += Vertices::Fw(sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(beta,a,sminputs,U)) * LoopFunctions::MFVV(pow(mnu[a],2), pow(sminputs.mW,2));
}
return a1r;
}
std::complex<double> A2L(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
std::complex<double> a2l = {0,0};
double mW = sminputs.mW;
for(int a=0; a<6; a++)
a2l += -2. * Vertices::Fw(sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(beta,a,sminputs,U)) * LoopFunctions::G1(pow(mnu[a],2), mW*mW, mW*mW) * ml[beta];
return a2l;
}
std::complex<double> A2R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
std::complex<double> a2r = {0,0};
double mW = sminputs.mW;
for(int a=0; a<6; a++)
a2r += -2. * Vertices::Fw(sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(beta,a,sminputs,U)) * LoopFunctions::G1(pow(mnu[a],2), mW*mW, mW*mW) * ml[alpha];
return a2r;
}
// Z penguins
std::complex<double> VZw2w4LL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
std::complex<double> vzll = {0,0};
for(int a=0; a<6; a++)
for(int c=0; c<3; c++)
{
// Use MZ for the renormalization scale Q
if(beta == c)
vzll += Vertices::EL(beta,c, sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),pow(sminputs.mW,2),sminputs.mZ))* pow(ml[alpha],2) / (pow(ml[alpha],2) - pow(ml[c],2));
if(alpha == c)
vzll += Vertices::EL(alpha,c, sminputs) * Vertices::VpL(beta,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),pow(sminputs.mW,2),sminputs.mZ))* pow(ml[beta],2) / (pow(ml[beta],2) - pow(ml[c],2));
}
return vzll;
}
std::complex<double> VZw2w4LR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return VZw2w4LL(alpha,beta,sminputs,U,ml,mnu);
}
std::complex<double> VZw2w4RR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
std::complex<double> vzrr = {0,0};
for(int a=0; a<6; a++)
for(int c=0; c<3; c++)
{
if(beta == c)
vzrr += Vertices::ER(beta,c,sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),pow(sminputs.mW,2),sminputs.mZ))* ml[c]*ml[alpha] / (pow(ml[alpha],2) - pow(ml[c],2));
if(alpha == c)
vzrr += Vertices::ER(alpha,c, sminputs) * Vertices::VpL(beta,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),pow(sminputs.mW,2),sminputs.mZ))* ml[c]*ml[beta] / (pow(ml[beta],2) - pow(ml[c],2));
}
return vzrr;
}
std::complex<double> VZw2w4RL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return VZw2w4RR(alpha, beta, sminputs, U, ml, mnu);
}
std::complex<double> VZw8LL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
std::complex<double> vzll = {0,0};
double mW = sminputs.mW;
// Use MZ as the renormalization scale Q
for(int a=0; a<6; a++)
vzll += Vertices::Zww(sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(beta,a,sminputs,U)) * (1. - 2.*(LoopFunctions::B0(mW*mW,mW*mW,sminputs.mZ) + 2.*LoopFunctions::C00(pow(mnu[a],2),mW*mW,mW*mW,sminputs.mZ) + LoopFunctions::C0(pow(mnu[a],2),mW*mW,mW*mW)*pow(mnu[a],2)));
return vzll;
}
std::complex<double> VZw8LR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
return VZw8LL(alpha, beta, sminputs, U, mnu);
}
std::complex<double> VZw10LL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
std::complex<double> vzll = {0,0};
double mW = sminputs.mW;
// Use MZ as the renormalization scale Q
for(int b=0; b<6; b++)
{
// Use different loop function in case that mnu[b] 0
if(mnu[b])
vzll += - Vertices::VpL(alpha,b,sminputs,U) * conj(Vertices::VpL(beta,b,sminputs,U)) * (2.* Vertices::VR(b,b,sminputs) * LoopFunctions::C0(pow(mnu[b],2),pow(mnu[b],2),mW*mW) * mnu[b] * mnu[b] + Vertices::VL(b,b,sminputs) * (1. - 2.*(LoopFunctions::B0(pow(mnu[b],2),pow(mnu[b],2),sminputs.mZ) - 2.*LoopFunctions::C00(pow(mnu[b],2),pow(mnu[b],2),mW*mW,sminputs.mZ) + LoopFunctions::C0(pow(mnu[b],2),pow(mnu[b],2),mW*mW)*mW*mW)));
else
vzll += - Vertices::VpL(alpha,b,sminputs,U) * conj(Vertices::VpL(beta,b,sminputs,U)) * Vertices::VL(b,b,sminputs) * (1. - 2.*(LoopFunctions::B02C00C0(pow(mnu[b],2),pow(mnu[b],2),mW*mW,sminputs.mZ)));
}
return vzll;
}
std::complex<double> VZw10LR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
return VZw10LL(alpha, beta, sminputs, U, mnu);
}
// Sum over Z penguins
std::complex<double> VZsumLL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return 1. / (16.*pow(pi,2)) * (VZw2w4LL(alpha, beta, sminputs, U, ml, mnu) + VZw8LL(alpha, beta, sminputs, U, mnu) + VZw10LL(alpha, beta, sminputs, U, mnu));
}
std::complex<double> VZsumLR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return 1. / (16.*pow(pi,2)) * (VZw2w4LR(alpha, beta, sminputs, U, ml, mnu) + VZw8LR(alpha, beta, sminputs, U, mnu) + VZw10LR(alpha, beta, sminputs, U, mnu));
}
std::complex<double> VZsumRL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return 1. / (16.*pow(pi,2)) * (VZw2w4RL(alpha, beta, sminputs, U, ml, mnu));
}
std::complex<double> VZsumRR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return 1. / (16.*pow(pi,2)) * (VZw2w4RR(alpha, beta, sminputs, U, ml, mnu));
}
// Scalar penguins
std::complex<double> Shw2w4LL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
std::complex<double> shll = {0,0};
double mW = sminputs.mW;
// Use mZ for the renormalisation scale Q
for(int a=0; a<6; a++)
for(int c=0; c<3; c++)
{
if(beta == c)
shll += - (Vertices::HL(beta,c,sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),mW*mW, sminputs.mZ)) * pow(ml[alpha],2))/(pow(ml[alpha],2) - pow(ml[c],2));
if(alpha == c)
shll += - (Vertices::HL(alpha,c,sminputs) * Vertices::VpL(beta,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),mW*mW, sminputs.mZ)) * pow(ml[beta],2))/(pow(ml[beta],2) - pow(ml[c],2));
}
return shll;
}
std::complex<double> Shw2w4LR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Shw2w4LL(alpha, beta, sminputs, U, ml, mnu);
}
std::complex<double> Shw2w4RR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
std::complex<double> shrr = {0,0};
double mW = sminputs.mW;
// Use mZ for the renormalisation scale Q
for(int a=0; a<6; a++)
for(int c=0; c<3; c++)
{
if(beta == c)
shrr += - (Vertices::HR(beta,c,sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),mW*mW, sminputs.mZ)) * ml[c]*ml[alpha])/(pow(ml[alpha],2) - pow(ml[c],2));
if(alpha == c)
shrr += - (Vertices::HR(alpha,c,sminputs) * Vertices::VpL(beta,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),mW*mW, sminputs.mZ)) * ml[c]*ml[beta])/(pow(ml[beta],2) - pow(ml[c],2));
}
return shrr;
}
std::complex<double> Shw2w4RL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Shw2w4RR(alpha, beta, sminputs, U, ml, mnu);
}
// Sum over scalar penguins
std::complex<double> ShsumLL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return 1. / (16.*pow(pi,2)) * Shw2w4LL(alpha, beta, sminputs, U, ml, mnu);
}
std::complex<double> ShsumLR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return 1. / (16.*pow(pi,2)) * Shw2w4LR(alpha, beta, sminputs, U, ml, mnu);
}
std::complex<double> ShsumRL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return 1. / (16.*pow(pi,2)) * Shw2w4RL(alpha, beta, sminputs, U, ml, mnu);
}
std::complex<double> ShsumRR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return 1. / (16.*pow(pi,2)) * Shw2w4RR(alpha, beta, sminputs, U, ml, mnu);
}
}
// Box contributions
namespace Boxes
{
std::complex<double> Vw4lLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
std::complex<double> vll = {0,0};
double mW = sminputs.mW;
for(int a=0; a<6; a++)
for(int c=0; c<6; c++)
vll += -4. * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(beta,a,sminputs,U)) * Vertices::VpL(gamma,c,sminputs,U) * conj(Vertices::VpL(delta,c,sminputs,U)) * (LoopFunctions::IC0D0(pow(mnu[c],2),mW*mW, mW*mW, pow(mnu[a],2)) - 3. * LoopFunctions::D27(pow(mnu[a],2),pow(mnu[c],2),mW*mW,mW*mW));
return vll;
}
std::complex<double> Vw8lLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
std::complex<double> vll = {0,0};
double mW = sminputs.mW;
for(int a=0; a<6; a++)
for(int c=0; c<6; c++)
vll += -2. * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(delta,c,sminputs,U)) * Vertices::VpL(gamma,a,sminputs,U) * conj(Vertices::VpL(beta,c,sminputs,U)) * mnu[a] * mnu[c] * LoopFunctions::D0(pow(mnu[a],2),pow(mnu[c],2),mW*mW,mW*mW);
return vll;
}
std::complex<double> Vw4lpLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
return Vw4lLL(alpha, delta, gamma, beta, sminputs, U, mnu);
}
std::complex<double> Vw8lpLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
return Vw8lLL(alpha, delta, gamma, beta, sminputs, U, mnu);
}
std::complex<double> Vw4dLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
std::complex<double> vll = {0,0};
double mW = sminputs.mW;
std::vector<double> mu = {sminputs.mU, sminputs.mCmC, sminputs.mT};
for(int a=0; a<6; a++)
for(int c=0; c<3; c++)
vll += -4.*Vertices::VpL(alpha,a,sminputs,U)*conj(Vertices::VpL(beta,a,sminputs,U))*Vertices::VuL(gamma,c,sminputs)*conj(Vertices::VuL(delta,c,sminputs))*(LoopFunctions::IC0D0(pow(mu[c],2),mW*mW, mW*mW, pow(mnu[a],2)) - 3.*LoopFunctions::D27(pow(mnu[a],2),pow(mu[c],2),mW*mW,mW*mW));
return vll;
}
std::complex<double> Vw4uLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
std::complex<double> vll = {0,0};
double mW = sminputs.mW;
std::vector<double> md = {sminputs.mD, sminputs.mS, sminputs.mBmB};
for(int a=0; a<6; a++)
for(int c=0; c<3; c++)
vll += 16.*Vertices::VpL(alpha,a,sminputs,U)*conj(Vertices::VpL(beta,a,sminputs,U))*Vertices::VuL(delta,c,sminputs)*conj(Vertices::VuL(gamma,c,sminputs))*LoopFunctions::D27(pow(mnu[a],2),pow(md[c],2),mW*mW,mW*mW);
return vll;
}
// Sum over boxes
std::complex<double> VsumlLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
return 1. / (16.*pow(pi,2)) *( Vw4lLL(alpha, beta, gamma, delta, sminputs, U, mnu) + Vw8lLL(alpha, beta, gamma, delta, sminputs, U, mnu) + Vw4lpLL(alpha, beta, gamma, delta, sminputs, U, mnu) + Vw8lpLL(alpha, beta, gamma, delta, sminputs, U, mnu));
}
std::complex<double> VsumdLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
return 1./(16.*pow(pi,2)) *Vw4dLL(alpha, beta, gamma, delta, sminputs, U, mnu);
}
std::complex<double> VsumuLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
return 1./(16.*pow(pi,2)) *Vw4uLL(alpha, beta, gamma, delta, sminputs, U, mnu);
}
} // Diagrams
// Form factors for LFV diagrams
namespace FormFactors
{
std::complex<double> K1R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> mnu)
{
double e = sqrt(4. * pi / sminputs.alphainv);
return 1. / (16*pow(pi,2)*e) * Penguins::A1R(alpha, beta, sminputs, U, mnu);
}
std::complex<double> K2L(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
double e = sqrt(4. * pi / sminputs.alphainv);
return 1. / (2. * 16.*pow(pi,2) * e * ml[alpha] ) * Penguins::A2L(alpha, beta, sminputs, U, ml, mnu);
}
std::complex<double> K2R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
double e = sqrt(4. * pi / sminputs.alphainv);
return 1. / (2. * 16.*pow(pi,2)* e * ml[alpha] ) * Penguins::A2R(alpha, beta, sminputs, U, ml, mnu);
}
std::complex<double> AVLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Penguins::VZsumLL(alpha,beta,sminputs,U,ml,mnu)*Vertices::EL(gamma,delta,sminputs) / pow(sminputs.mZ,2) + Boxes::VsumlLL(alpha,beta,gamma,delta,sminputs,U,mnu);
}
std::complex<double> AVLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Penguins::VZsumLR(alpha,beta,sminputs,U,ml,mnu)*Vertices::ER(gamma,delta,sminputs) / pow(sminputs.mZ,2);
}
std::complex<double> AVRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Penguins::VZsumRL(alpha,beta,sminputs,U,ml,mnu)*Vertices::EL(gamma,delta,sminputs) / pow(sminputs.mZ,2);
}
std::complex<double> AVRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Penguins::VZsumRR(alpha,beta,sminputs,U,ml,mnu)*Vertices::ER(gamma,delta,sminputs) / pow(sminputs.mZ,2);
}
std::complex<double> ASLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu, double mh)
{
return Penguins::ShsumLL(alpha,beta,sminputs,U,ml,mnu)*Vertices::HL(gamma,delta,sminputs) / pow(mh,2);
}
std::complex<double> ASLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu, double mh)
{
return Penguins::ShsumLR(alpha,beta,sminputs,U,ml,mnu)*Vertices::HR(gamma,delta,sminputs) / pow(mh,2);
}
std::complex<double> ASRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu, double mh)
{
return Penguins::ShsumRL(alpha,beta,sminputs,U,ml,mnu)*Vertices::HL(gamma,delta,sminputs) / pow(mh,2);
}
std::complex<double> ASRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu, double mh)
{
return Penguins::ShsumRR(alpha,beta,sminputs,U,ml,mnu)*Vertices::HR(gamma,delta,sminputs) / pow(mh,2);
}
std::complex<double> BVLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Penguins::VZsumLL(alpha,beta,sminputs,U,ml,mnu)*Vertices::DL(gamma,delta,sminputs) / pow(sminputs.mZ,2) + Boxes::VsumdLL(alpha,beta,gamma,delta,sminputs,U,mnu);
}
std::complex<double> BVLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Penguins::VZsumLR(alpha,beta,sminputs,U,ml,mnu)*Vertices::DR(gamma,delta,sminputs) / pow(sminputs.mZ,2);
}
std::complex<double> BVRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Penguins::VZsumRL(alpha,beta,sminputs,U,ml,mnu)*Vertices::DL(gamma,delta,sminputs) / pow(sminputs.mZ,2);
}
std::complex<double> BVRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Penguins::VZsumRR(alpha,beta,sminputs,U,ml,mnu)*Vertices::DR(gamma,delta,sminputs) / pow(sminputs.mZ,2);
}
std::complex<double> BSLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu, double mh)
{
return Penguins::ShsumLL(alpha,beta,sminputs,U,ml,mnu)*Vertices::HdL(gamma,delta,sminputs) / pow(mh,2);
}
std::complex<double> BSLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu, double mh)
{
return Penguins::ShsumLR(alpha,beta,sminputs,U,ml,mnu)*Vertices::HdR(gamma,delta,sminputs) / pow(mh,2);
}
std::complex<double> BSRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu, double mh)
{
return Penguins::ShsumRL(alpha,beta,sminputs,U,ml,mnu)*Vertices::HdL(gamma,delta,sminputs) / pow(mh,2);
}
std::complex<double> BSRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu, double mh)
{
return Penguins::ShsumRR(alpha,beta,sminputs,U,ml,mnu)*Vertices::HdR(gamma,delta,sminputs) / pow(mh,2);
}
std::complex<double> CVLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Penguins::VZsumLL(alpha,beta,sminputs,U,ml,mnu)*Vertices::UL(gamma,delta,sminputs) / pow(sminputs.mZ,2) + Boxes::VsumuLL(alpha,beta,gamma,delta,sminputs,U,mnu);
}
std::complex<double> CVLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Penguins::VZsumLR(alpha,beta,sminputs,U,ml,mnu)*Vertices::UR(gamma,delta,sminputs) / pow(sminputs.mZ,2);
}
std::complex<double> CVRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Penguins::VZsumRL(alpha,beta,sminputs,U,ml,mnu)*Vertices::UL(gamma,delta,sminputs) / pow(sminputs.mZ,2);
}
std::complex<double> CVRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu)
{
return Penguins::VZsumRR(alpha,beta,sminputs,U,ml,mnu)*Vertices::UR(gamma,delta,sminputs) / pow(sminputs.mZ,2);
}
std::complex<double> CSLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu, double mh)
{
return Penguins::ShsumLL(alpha,beta,sminputs,U,ml,mnu)*Vertices::HuL(gamma,delta,sminputs) / pow(mh,2);
}
std::complex<double> CSLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu, double mh)
{
return Penguins::ShsumLR(alpha,beta,sminputs,U,ml,mnu)*Vertices::HuR(gamma,delta,sminputs) / pow(mh,2);
}
std::complex<double> CSRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu, double mh)
{
return Penguins::ShsumRL(alpha,beta,sminputs,U,ml,mnu)*Vertices::HuL(gamma,delta,sminputs) / pow(mh,2);
}
std::complex<double> CSRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<std::complex<double>,3,6> U, std::vector<double> ml, std::vector<double> mnu, double mh)
{
return Penguins::ShsumRR(alpha,beta,sminputs,U,ml,mnu)*Vertices::HuR(gamma,delta,sminputs) / pow(mh,2);
}
} // Form Factors
}
}
#endif //#defined __flav_loop_functions_hpp__
Updated on 2023-06-26 at 21:36:54 +0000