2double CLASSNAME::get_partial_width<Higgs, bar<UpTypeQuark>::type, UpTypeQuark>(
3 const context_base& context,
4 typename field_indices<Higgs>::type
const& indexIn,
5 typename field_indices<UpTypeQuark>::type
const& indexOut1,
6 typename field_indices<UpTypeQuark>::type
const& indexOut2)
11 const auto indices =
concatenate(indexOut1, indexOut2, indexIn);
12 const auto HQQbarVertexDR = Vertex<bar<UpTypeQuark>::type, UpTypeQuark, Higgs>::evaluate(indices, context);
14 const double mHOS = context.physical_mass<Higgs>(indexIn);
15 const double flux = 1./(2.*mHOS);
17 static constexpr double color_factor = squared_color_generator<Higgs, bar<UpTypeQuark>::type, UpTypeQuark>();
19 if(!boost::range::equal(indexOut1, indexOut2)) {
20 if (!
is_zero(HQQbarVertexDR.left()) || !
is_zero(HQQbarVertexDR.right())) {
21 const double muqOS1 = context.physical_mass<UpTypeQuark>(indexOut1);
22 const double muqOS2 = context.physical_mass<UpTypeQuark>(indexOut2);
23 const auto xOS1 =
Sqr(muqOS1/mHOS);
24 const auto xOS2 =
Sqr(muqOS2/mHOS);
25 const double phase_space = 1./(8.*
Pi) * std::sqrt(
KallenLambda(1., xOS1, xOS2));
26 return flux * phase_space * color_factor * amplitude_squared<Higgs, bar<UpTypeQuark>::type, UpTypeQuark>(context, indexIn, indexOut1, indexOut2);
31 const double muqDR = context.mass<UpTypeQuark>(indexOut1);
32 const double muqOS = context.physical_mass<UpTypeQuark>(indexOut1);
34 throw std::runtime_error(
36 +
": Up-type quark cannot be massless. Aborting."
39 const auto xOS =
Sqr(muqOS/mHOS);
40 const auto xDR =
Sqr(muqDR/mHOS);
42 const auto HQQbarVertexDR_S = 0.5*(HQQbarVertexDR.left() + HQQbarVertexDR.right());
43 const auto HQQbarVertexDR_P = 0.5*(HQQbarVertexDR.right() - HQQbarVertexDR.left());
48 if (4.*std::max(xDR, xOS) > 1.) {
51 const auto betaOS = std::sqrt(1.-4.*xOS);
52 const auto betaDR = std::sqrt(1.-4.*xDR);
54 const double phase_spaceDR = 1./(8.*
Pi) * std::sqrt(
KallenLambda(1., xDR, xDR));
55 const double phase_spaceOS = 1./(8.*
Pi) * std::sqrt(
KallenLambda(1., xOS, xOS));
57 double amp2DR_S =
Sqr(mHOS) *
Sqr(betaDR) *
59 double amp2OS_S =
Sqr(mHOS) *
Sqr(betaOS) *
64 if (info::is_CP_violating_Higgs_sector) {
65 amp2DR_P =
Sqr(mHOS) *
67 amp2OS_P =
Sqr(mHOS) *
75 +
": Cannot determine the number of active quark flavours. Disabling higher-order corrections."
79 if (
static_cast<int>(flexibledecay_settings.get(FlexibleDecay_settings::include_higher_order_corrections)) > 0 && Nf >= 5) {
81 double Y_conversion = 1.;
84 auto qedqcd_ = qedqcd;
87 Y_conversion =
Sqr(qedqcd_.displayUpQuarkRunningMass(indexOut1.at(0))/muqDR);
91 alpha_s_red = get_alphas(context)/
Pi;
94 throw std::runtime_error(
96 +
": Cannot determine the number of active quark flavours"
99 double deltaqq_QCD_DR_S =
calc_Deltaqq(alpha_s_red, Nf, flexibledecay_settings);
100 double deltaqq_QCD_DR_P = deltaqq_QCD_DR_S;
103 const double alpha_red = get_alpha(context)/
Pi;
104 const double deltaqq_QED_DR = 17./4.*
Sqr(UpTypeQuark::electricCharge)*alpha_red;
107 2.*(1. - 10.*xDR)/(1-4.*xDR)*(4./3. -
std::log(xDR))*alpha_s_red +
110 const double deltaqq_QCD_OS_S =
113 const double deltaqq_QED_OS_S =
114 alpha_red *
Sqr(UpTypeQuark::electricCharge) *
calc_DeltaH(betaOS);
116 double deltaqq_QCD_OS_P = 0.;
117 double deltaqq_QED_OS_P = 0.;
119 if (info::is_CP_violating_Higgs_sector) {
121 2.*(1. - 6.*xDR)/(1-4.*xDR)*(4./3. -
std::log(xDR))*alpha_s_red +
131 double deltaPhi2_S = 0.;
132 double deltaPhi2_P = 0.;
133 double deltaqq_QCDxQED_DR = 0.;
134 if (
static_cast<int>(flexibledecay_settings.get(FlexibleDecay_settings::include_higher_order_corrections)) > 1) {
137 if ((indexOut1.at(0) < 2 || indexOut2.at(0) < 2)) {
138 const double mtpole = qedqcd.displayPoleMt();
142 const auto Httindices =
concatenate(std::array<int, 1> {2}, std::array<int, 1> {2}, indexIn);
143 const auto Httbar = Vertex<bar<UpTypeQuark>::type, UpTypeQuark, Higgs>::evaluate(Httindices, context);
144 const auto gbHoVEV = HQQbarVertexDR_S/muqDR;
146 const auto Httbar_S = 0.5*(Httbar.left() + Httbar.right());
147 const auto gtHoVEV = Httbar_S/context.mass<UpTypeQuark>({2});
148 deltaPhi2_S =
Sqr(alpha_s_red) * std::real(gtHoVEV/gbHoVEV) * (8/3. -
Sqr(
Pi/3.) - 2.0/3.0*lt + 1.0/9.0*
Sqr(lq));
150 if (info::is_CP_violating_Higgs_sector) {
151 const auto CSuu = HQQbarVertexDR_S/muqDR;
153 const auto Httbar_P = 0.5*(Httbar.right() - Httbar.left());
154 const auto CStu = Httbar_P/context.mass<Fu>({2});
155 deltaPhi2_P =
Sqr(alpha_s_red) * std::real(CStu/CSuu) * (23/6. - lt + 1.0/6.0*
Sqr(lq));
161 amp2DR_S *= Y_conversion*(1. + deltaqq_QCD_DR_S + deltaqq_QED_DR + deltaqq_QCDxQED_DR + deltaPhi2_S);
162 amp2DR_P *= Y_conversion*(1. + deltaqq_QCD_DR_P + deltaqq_QED_DR + deltaqq_QCDxQED_DR + deltaPhi2_P);
163 amp2OS_S *= 1. + deltaqq_QCD_OS_S + deltaqq_QED_OS_S;
164 amp2OS_P *= 1. + deltaqq_QCD_OS_P + deltaqq_QED_OS_P;
168 const double result_DR_S =
169 flux * color_factor * phase_spaceDR * amp2DR_S;
170 const double result_DR_P =
171 flux * color_factor * phase_spaceDR * amp2DR_P;
173 const double result_OS_S =
174 flux * color_factor * phase_spaceOS * amp2OS_S;
175 const double result_OS_P =
176 flux * color_factor * phase_spaceOS * amp2OS_P;
178 result_S = (1-4.*xOS)*result_DR_S + 4*xOS*result_OS_S;
179 result_P = (1-4.*xOS)*result_DR_P + 4*xOS*result_OS_P;
182 return result_S + result_P;
bool is_zero(double x) noexcept
constexpr T norm(const Complex< T > &z) noexcept
static constexpr double deltaqq_QCDxQED
double calc_DeltaAh(double b) noexcept
std::string create_process_string(std::array< int, FieldIn::numberOfFieldIndices > const in, std::array< int, FieldOut1::numberOfFieldIndices > const out1, std::array< int, FieldOut2::numberOfFieldIndices > const out2)
unsigned int number_of_active_flavours(softsusy::QedQcd const &qedqcd, double m) noexcept
double calc_DeltaH(double b) noexcept
Eq.(2.6) of hep-ph/0503173.
double calc_Deltaqq(double alpha_s_red, double Nf, FlexibleDecay_settings const &settings) noexcept
Eq.(2.11) of hep-ph/0503173, 2-loop and higher order.
constexpr std::complex< T > Sqr(const std::complex< T > &a) noexcept
detail::result_of::concatenate< Args... >::type concatenate(Args &&... args)
T KallenLambda(T x, T y, T z) noexcept
Complex< T > log(const Complex< T > &z) noexcept
static constexpr double Pi