2double CLASSNAME::get_partial_width<PseudoscalarHiggs, bar<UpTypeQuark>::type, UpTypeQuark>(
3 const context_base& context,
4 typename field_indices<PseudoscalarHiggs>::type
const& indexIn,
5 typename field_indices<UpTypeQuark>::type
const& indexOut1,
6 typename field_indices<UpTypeQuark>::type
const& indexOut2)
8 if (indexIn.at(0) < info::number_of_neutral_goldstones) {
9 throw OutOfBoundsError(
"Error in " +
create_process_string<PseudoscalarHiggs,bar<UpTypeQuark>::type, UpTypeQuark>(indexIn, indexOut1, indexOut2) +
" decay. Decaying particle is a Goldstone.");
15 const auto indices =
concatenate(indexOut1, indexOut2, indexIn);
16 const auto AhQQbarVertexDR = Vertex<bar<UpTypeQuark>::type, UpTypeQuark, PseudoscalarHiggs>::evaluate(indices, context);
18 const double mAhOS = context.physical_mass<PseudoscalarHiggs>(indexIn);
19 const double flux = 1./(2.*mAhOS);
21 static constexpr double color_factor = squared_color_generator<PseudoscalarHiggs, bar<UpTypeQuark>::type, UpTypeQuark>();
24 if(!boost::range::equal(indexOut1, indexOut2)) {
25 if (!
is_zero(AhQQbarVertexDR.left()) || !
is_zero(AhQQbarVertexDR.right())) {
26 const double muqOS1 = context.physical_mass<UpTypeQuark>(indexOut1);
27 const double muqOS2 = context.physical_mass<UpTypeQuark>(indexOut2);
28 const auto xOS1 =
Sqr(muqOS1/mAhOS);
29 const auto xOS2 =
Sqr(muqOS2/mAhOS);
30 const double phase_space = 1./(8.*
Pi) * std::sqrt(
KallenLambda(1., xOS1, xOS2));
31 return flux * phase_space * color_factor * amplitude_squared<PseudoscalarHiggs, bar<UpTypeQuark>::type, UpTypeQuark>(context, indexIn, indexOut1, indexOut2);
36 const double muqDR = context.mass<UpTypeQuark>(indexOut1);
37 const double muqOS = context.physical_mass<UpTypeQuark>(indexOut1);
39 throw std::runtime_error(
40 create_process_string<PseudoscalarHiggs,bar<UpTypeQuark>::type, UpTypeQuark>(indexIn, indexOut1, indexOut2)
41 +
": Up-type quark cannot be massless. Aborting."
44 const auto xOS =
Sqr(muqOS/mAhOS);
45 const auto xDR =
Sqr(muqDR/mAhOS);
49 const auto AhQQbarVertexDR_P = 0.5*(AhQQbarVertexDR.right() - AhQQbarVertexDR.left());
52 if (4.*std::max(xDR, xOS) > 1.) {
55 const auto betaOS = std::sqrt(1.-4.*xOS);
56 const auto betaDR = std::sqrt(1.-4.*xDR);
58 const double phase_spaceDR = 1./(8.*
Pi) * std::sqrt(
KallenLambda(1., xDR, xDR));
59 const double phase_spaceOS = 1./(8.*
Pi) * std::sqrt(
KallenLambda(1., xOS, xOS));
63 amp2DR_P =
Sqr(mAhOS) *
65 amp2OS_P =
Sqr(mAhOS) *
71 create_process_string<PseudoscalarHiggs,bar<UpTypeQuark>::type, UpTypeQuark>(indexIn, indexOut1, indexOut2)
72 +
": Cannot determine the number of active quark flavours. Disabling higher-order corrections."
76 if (
static_cast<int>(flexibledecay_settings.get(FlexibleDecay_settings::include_higher_order_corrections)) > 0 && Nf >= 5) {
78 double Y_conversion = 1.;
81 auto qedqcd_ = qedqcd;
84 Y_conversion =
Sqr(qedqcd_.displayUpQuarkRunningMass(indexOut1.at(0))/muqDR);
88 alpha_s_red = get_alphas(context)/
Pi;
91 throw std::runtime_error(
92 create_process_string<PseudoscalarHiggs,bar<UpTypeQuark>::type, UpTypeQuark>(indexIn, indexOut1, indexOut2)
93 +
": Cannot determine the number of active quark flavours"
96 double deltaqq_QCD_DR_P =
98 + 2.*(1. - 6.*xDR)/(1-4.*xDR)*(4./3. -
std::log(xDR))*alpha_s_red
102 const double alpha_red = get_alpha(context)/
Pi;
103 const double deltaqq_QED_DR = 17./4.*
Sqr(UpTypeQuark::electricCharge)*alpha_red;
105 const double deltaqq_QCD_OS_P =
108 const double deltaqq_QED_OS_P =
111 double deltaqq_QCDxQED_DR = 0.;
112 double deltaPhi2_P = 0.;
113 if (
static_cast<int>(flexibledecay_settings.get(FlexibleDecay_settings::include_higher_order_corrections)) > 1) {
116 if ((indexOut1.at(0) < 2 || indexOut2.at(0) < 2)) {
117 const double mtpole = qedqcd.displayPoleMt();
121 const auto Ahttindices =
concatenate(std::array<int, 1> {2}, std::array<int, 1> {2}, indexIn);
122 const auto Ahttbar = Vertex<bar<UpTypeQuark>::type, UpTypeQuark, PseudoscalarHiggs>::evaluate(Ahttindices, context);
123 const auto CSuu = AhQQbarVertexDR_P/context.mass<UpTypeQuark>(indexOut1);
125 const auto Ahttbar_P = 0.5*(Ahttbar.right() - Ahttbar.left());
126 const auto CStu = Ahttbar_P/context.mass<Fu>({2});
127 deltaPhi2_P =
Sqr(alpha_s_red) * std::real(CStu/CSuu) * (23/6. - lt + 1.0/6.0*
Sqr(lq));
132 amp2DR_P *= Y_conversion*(1. + deltaqq_QCD_DR_P + deltaqq_QED_DR + deltaqq_QCDxQED_DR + deltaPhi2_P);
133 amp2OS_P *= 1. + deltaqq_QCD_OS_P + deltaqq_QED_OS_P;
138 flux * color_factor * phase_spaceDR * amp2DR_P;
141 flux * color_factor * phase_spaceOS * amp2OS_P;
143 result = (1-4.*xOS)*result_DR + 4*xOS*result_OS;
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_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