Hyperiso 1.0.3
Modular flavour-physics calculations, Wilson coefficients and statistical inference
Loading...
Searching...
No Matches
NMSSMScalarMatching.cpp
Go to the documentation of this file.
1#include "NMSSMScalarMatching.h"
2
3#include <array>
4#include <cmath>
5#include <complex>
6#include <sstream>
7#include <stdexcept>
8#include <string>
9
10#include "Logger.h"
11#include "special_SUSY.h"
12#include "scalar.h"
13#include "ParameterProvider.h"
14
15namespace nmssm_scalar_matching {
16namespace {
17
18constexpr double SQRT2 = 1.4142135623730950488;
19constexpr double TINY = 1.0e-14;
20
21using M22 = std::array<std::array<double, 2>, 2>;
22using M23 = std::array<std::array<double, 3>, 2>;
23using M33 = std::array<std::array<double, 3>, 3>;
24using A222 = std::array<std::array<std::array<double, 2>, 2>, 2>;
25using A322 = std::array<std::array<std::array<double, 2>, 2>, 3>;
26using A333 = std::array<std::array<std::array<double, 3>, 3>, 3>;
27using A3322 = std::array<std::array<std::array<std::array<double, 2>, 2>, 3>, 3>;
28
29inline double real_value(const scalar_t& value) {
30 return static_cast<complex_t>(value).real();
31}
32
33bool exists(ParameterType type, const std::string& block, const LhaID& id) {
34 return ParameterProvider(type).exists(block, id);
35}
36
37double get(ParameterType type, const std::string& block, const LhaID& id) {
38 return real_value(ParameterProvider(type)(block, id));
39}
40
41double get_required(ParameterType type, const std::string& block, const LhaID& id,
42 const std::string& label) {
43 if (!exists(type, block, id)) {
44 std::ostringstream os;
45 os << "NMSSM scalar matching: missing " << label << " (" << block << ", " << id.to_string() << ")";
46 throw std::runtime_error(os.str());
47 }
48 return get(type, block, id);
49}
50
51double get_nmssm_parameter(int nmssmrun_id, int extpar_id, const std::string& label,
52 bool allow_hmix_mu = false) {
53 if (exists(ParameterType::BSM, "NMSSMRUN", LhaID(nmssmrun_id))) {
54 return get(ParameterType::BSM, "NMSSMRUN", LhaID(nmssmrun_id));
55 }
56 if (exists(ParameterType::PASSTHROUGH, "EXTPAR", LhaID(extpar_id))) {
57 return get(ParameterType::PASSTHROUGH, "EXTPAR", LhaID(extpar_id));
58 }
59 if (allow_hmix_mu && exists(ParameterType::BSM, "HMIX", LhaID(1))) {
60 return get(ParameterType::BSM, "HMIX", LhaID(1));
61 }
62 std::ostringstream os;
63 os << "NMSSM scalar matching: missing " << label
64 << " (NMSSMRUN " << nmssmrun_id << " or EXTPAR " << extpar_id << ")";
65 throw std::runtime_error(os.str());
66}
67
68inline double sqr(double x) { return x * x; }
69inline double delta(int i, int j) { return i == j ? 1.0 : 0.0; }
70
71void require_nonzero(double value, const std::string& label) {
72 if (std::abs(value) < TINY) {
73 throw std::runtime_error("NMSSM scalar matching: zero/near-zero " + label);
74 }
75}
76
77} // namespace
78
79bool is_active() {
80 const ParameterProvider bsm(ParameterType::BSM);
81 return (bsm.exists("MASS", LhaID(45)) && std::abs(real_value(bsm("MASS", LhaID(45)))) > TINY)
82 || (bsm.exists("MASS", LhaID(46)) && std::abs(real_value(bsm("MASS", LhaID(46)))) > TINY);
83}
84
85Result compute(const ParamSrc& src, int lepton_mass_slot) {
86 if (!is_active()) {
87 return {};
88 }
89
90 const double lambda = get_nmssm_parameter(1, 61, "lambda");
91 const double kappa_nmssm = get_nmssm_parameter(2, 62, "kappa");
92 const double A_lambda = get_nmssm_parameter(3, 63, "A_lambda");
93 const double lambda_s = get_nmssm_parameter(5, 65, "lambda*<S> = mu_eff", true);
94 require_nonzero(lambda, "lambda");
95
96 const double singlet_vev = lambda_s / lambda;
97 require_nonzero(singlet_vev, "<S>");
98
99 const double gf = get_required(ParameterType::SM, "SMINPUTS", LhaID(2), "G_F");
100 const double v = std::sqrt(1.0 / std::sqrt(2.0) / gf);
101 const double vdeltam_den = SQRT2 * A_lambda + kappa_nmssm * singlet_vev;
102 require_nonzero(vdeltam_den, "sqrt(2) A_lambda + kappa <S>");
103 const double v_deltam_s = v / singlet_vev
104 * (SQRT2 * A_lambda - 2.0 * kappa_nmssm * singlet_vev) / vdeltam_den;
105
106 const double mW = real_value(src.get_val(ParameterType::SM, "MASS", 24));
107 const double mZ = get_required(ParameterType::SM, "SMINPUTS", LhaID(4), "M_Z");
108 const double mHc = real_value(src.get_val(ParameterType::BSM, "MASS", 37));
109 const double sw2 = real_value(src.get_val(ParameterType::WILSON, "WPARAM_SI_SM", 4));
110 const double mb_muW = real_value(src.get_val(ParameterType::WILSON, "WPARAM_MATCH_SM", {5, 1}));
111 const double mt_muW = real_value(src.get_val(ParameterType::WILSON, "WPARAM_MATCH_SM", 6));
112 const double g2 = real_value(src.get_val(ParameterType::SM, "GAUGE", 2));
113 const double tanb = real_value(src.get_val(ParameterType::BSM, "MINPAR", 3));
114 const double epsfac = real_value(src.get_val(ParameterType::WILSON, "EPSILON_SUSY", 5));
115 const double ml = real_value(src.get_val(ParameterType::WILSON, "WPARAM_SI_SM", lepton_mass_slot));
116 const double Qmatch = real_value(src.get_val(ParameterType::WILSON, "EW_SCALE", 1));
117
118 require_nonzero(mW, "M_W");
119 require_nonzero(mHc, "M_H+");
120 require_nonzero(g2, "g2");
121 require_nonzero(epsfac, "epsilon factor");
122
123 const std::array<double, 3> mh = {
124 get_required(ParameterType::SM, "MASS", LhaID(25), "MASS 25 (h1)"),
125 get_required(ParameterType::BSM, "MASS", LhaID(35), "MASS 35 (h2)"),
126 get_required(ParameterType::BSM, "MASS", LhaID(45), "MASS 45 (h3)")
127 };
128 const std::array<double, 2> ma = {
129 get_required(ParameterType::BSM, "MASS", LhaID(36), "MASS 36 (a1)"),
130 get_required(ParameterType::BSM, "MASS", LhaID(46), "MASS 46 (a2)")
131 };
132 const std::array<double, 3> mstop = {
133 get_required(ParameterType::BSM, "MASS", LhaID(2000002), "MASS 2000002 (u_R-like squark)"),
134 get_required(ParameterType::BSM, "MASS", LhaID(1000006), "MASS 1000006 (stop 1)"),
135 get_required(ParameterType::BSM, "MASS", LhaID(2000006), "MASS 2000006 (stop 2)")
136 };
137 const double m_snutau = get_required(ParameterType::BSM, "MASS", LhaID(1000016), "MASS 1000016 (tau sneutrino)");
138 const double Au = get_required(ParameterType::BSM, "AU", LhaID(1, 1), "AU(1,1)");
139
140 std::array<double, 2> mch {};
141 for (int j = 0; j < 2; ++j) {
142 mch[j] = real_value(src.get_val(ParameterType::WILSON, "WPARAM_SI_BSM", {13, j}));
143 require_nonzero(mch[j], "chargino mass");
144 }
145
146 M22 U {}, V {}, stopmix {};
147 M33 H {}, TU {};
148 M23 A {};
149 for (int i = 0; i < 2; ++i) {
150 for (int j = 0; j < 2; ++j) {
151 U[i][j] = get_required(ParameterType::BSM, "UMIX", LhaID(i + 1, j + 1), "UMIX");
152 V[i][j] = get_required(ParameterType::BSM, "VMIX", LhaID(i + 1, j + 1), "VMIX");
153 stopmix[i][j] = get_required(ParameterType::BSM, "STOPMIX", LhaID(i + 1, j + 1), "STOPMIX");
154 }
155 }
156 for (int a = 0; a < 3; ++a) {
157 for (int c = 0; c < 3; ++c) {
158 H[a][c] = get_required(ParameterType::BSM, "NMHMIX", LhaID(a + 1, c + 1), "NMHMIX");
159 }
160 }
161 for (int a = 0; a < 2; ++a) {
162 for (int c = 0; c < 3; ++c) {
163 A[a][c] = get_required(ParameterType::BSM, "NMAMIX", LhaID(a + 1, c + 1), "NMAMIX");
164 }
165 }
166
167 TU[0][0] = 1.0;
168 for (int i = 0; i < 2; ++i) {
169 for (int j = 0; j < 2; ++j) {
170 TU[i + 1][j + 1] = stopmix[i][j];
171 }
172 }
173
174 const double beta = std::atan(tanb);
175 const double sinb = std::sin(beta);
176 const double cosb = std::cos(beta);
177 require_nonzero(sinb, "sin(beta)");
178 require_nonzero(cosb, "cos(beta)");
179 const double vu = std::sqrt(sqr(sinb) / std::sqrt(2.0) / gf);
180 const double vd = vu / tanb;
181
182 A222 R {};
183 A322 Q {};
184 A3322 G1 {};
185 A333 T2 {};
186 std::array<std::array<std::array<double, 3>, 3>, 2> T1 {};
187
188 for (int a = 0; a < 2; ++a) {
189 for (int l = 0; l < 2; ++l) {
190 for (int j = 0; j < 2; ++j) {
191 R[a][l][j] = -g2 / SQRT2
192 * (A[a][0] * U[1][l] * V[1][j] + A[a][1] * U[0][l] * V[1][j])
193 - lambda / SQRT2 * A[a][2] * U[1][l] * V[1][j];
194 }
195 }
196 }
197 for (int a = 0; a < 3; ++a) {
198 for (int l = 0; l < 2; ++l) {
199 for (int j = 0; j < 2; ++j) {
200 Q[a][l][j] = g2 / SQRT2
201 * (H[a][0] * U[1][l] * V[1][j] + H[a][1] * U[0][l] * V[1][j])
202 - lambda / SQRT2 * H[a][2] * U[1][l] * V[1][j];
203 }
204 }
205 }
206 for (int i = 0; i < 3; ++i) {
207 for (int k = 0; k < 3; ++k) {
208 for (int j = 0; j < 2; ++j) {
209 for (int l = 0; l < 2; ++l) {
210 G1[i][k][j][l] = (TU[i][1] * TU[k][1] - delta(i, 0) * delta(k, 0)) * V[0][l] * U[1][j]
211 - mt_muW / SQRT2 / sinb / mW * TU[i][2] * TU[k][1] * V[1][l] * U[1][j];
212 }
213 }
214 }
215 }
216 for (int a = 0; a < 2; ++a) {
217 for (int i = 0; i < 3; ++i) {
218 for (int k = 0; k < 3; ++k) {
219 T1[a][i][k] = (TU[i][2] * TU[k][1] - TU[i][1] * TU[k][2])
220 * ((lambda / SQRT2 * (vd * A[a][2] + singlet_vev * A[a][0])) - Au * A[a][1]);
221 }
222 }
223 }
224 for (int a = 0; a < 3; ++a) {
225 for (int i = 0; i < 3; ++i) {
226 for (int k = 0; k < 3; ++k) {
227 T2[a][i][k] = -mt_muW / (2.0 * mW)
228 * (2.0 * mt_muW * H[a][1] * (TU[i][1] * TU[k][1] + TU[i][2] * TU[k][2])
229 + ((lambda / SQRT2 * (vd * H[a][2] + singlet_vev * H[a][0])) + Au * H[a][1])
230 * (TU[i][2] * TU[k][1] + TU[i][1] * TU[k][2]))
231 + mZ / 2.0 / std::sqrt(1.0 - sw2) * (1.0 - 4.0 / 3.0 * sw2) * H[a][1]
232 * (TU[i][0] * TU[k][0] + TU[i][1] * TU[k][1])
233 + 2.0 / 3.0 * mW * sw2 / (1.0 - sw2) * H[a][1] * TU[i][2] * TU[k][2];
234 }
235 }
236 }
237
238 complex_t cq1_h = 0.0;
239 for (int a = 0; a < 3; ++a) {
240 require_nonzero(mh[a], "CP-even Higgs mass");
241 cq1_h += (sqr(mHc / mW) * sqr(H[a][0]) * f30(sqr(mHc / mt_muW), sqr(mW / mt_muW))
242 + sqr(mt_muW) * sqr(mh[a]) / sqr(mW) / sqr(mHc)
243 * f30(sqr(mt_muW / mHc), sqr(mt_muW / mW))) / sqr(mh[a]);
244 }
245 cq1_h *= -ml / 4.0 * sqr(tanb);
246
247 complex_t cq2_h = 0.0;
248 for (int a = 0; a < 2; ++a) {
249 require_nonzero(ma[a], "CP-odd Higgs mass");
250 cq2_h += ((sqr(mHc / mW) * sqr(A[a][0]) + delta(a, 1) * A[a][0])
251 * f30(sqr(mHc / mt_muW), sqr(mW / mt_muW))
252 + sqr(mt_muW) * sqr(ma[a]) / sqr(mW) / sqr(mHc)
253 * f30(sqr(mt_muW / mHc), sqr(mt_muW / mW))) / sqr(ma[a]);
254 }
255 cq2_h *= ml / 4.0 * sqr(tanb);
256
257 const complex_t I(0.0, 1.0);
258 const complex_t ca_h = -I * lambda * A_lambda / g2 / mW * tanb
259 * f30(sqr(mHc / mt_muW), sqr(mW / mt_muW));
260
261 complex_t cq1_c = 0.0;
262 for (int a = 0; a < 3; ++a) {
263 for (int i = 0; i < 3; ++i) {
264 for (int k = 0; k < 3; ++k) {
265 for (int j = 0; j < 2; ++j) {
266 for (int l = 0; l < 2; ++l) {
267 cq1_c += G1[i][k][j][l] / sqr(mh[a]) * (
268 SQRT2 * sqr(H[a][0]) * mch[j] / mW / cosb * delta(i, k) * delta(l, j)
269 * f80(sqr(mstop[i] / mch[j]))
270 - 2.0 * SQRT2 * H[a][0] / g2 * delta(i, k)
271 * (Q[a][l][j] * f40(sqr(mstop[i] / mch[l]), sqr(mch[j] / mch[l]))
272 + mch[j] / mch[l] * Q[a][j][l]
273 * f30(sqr(mstop[i] / mch[l]), sqr(mch[j] / mch[l])))
274 + 2.0 * SQRT2 * H[a][0] * T2[a][i][k] * mch[j] / sqr(mstop[k]) * delta(l, j)
275 * f30(sqr(mstop[i] / mstop[k]), sqr(mch[j] / mstop[k]))
276 + sqr(mh[a] / mch[j]) * delta(i, k)
277 * (U[1][j] * V[0][l]
278 * f50(sqr(mstop[i] / mch[j]), sqr(mch[l] / mch[j]), sqr(m_snutau / mch[l]))
279 - mch[l] / mch[j] * U[1][l] * V[0][j]
280 * f60(sqr(mstop[i] / mch[j]), sqr(mch[l] / mch[j]), sqr(m_snutau / mch[l])))
281 );
282 }
283 }
284 }
285 }
286 }
287 cq1_c *= ml / 4.0 * sqr(tanb);
288
289 complex_t cq2_c = 0.0;
290 for (int a = 0; a < 2; ++a) {
291 for (int i = 0; i < 3; ++i) {
292 for (int k = 0; k < 3; ++k) {
293 for (int j = 0; j < 2; ++j) {
294 for (int l = 0; l < 2; ++l) {
295 cq2_c += G1[i][k][j][l] / sqr(ma[a]) * (
296 SQRT2 * sqr(A[a][0]) * mch[j] / mW / cosb * delta(i, k) * delta(l, j)
297 * f80(sqr(mstop[i] / mch[j]))
298 - 2.0 * SQRT2 * A[a][0] / g2 * delta(i, k)
299 * (-R[a][l][j] * f40(sqr(mstop[i] / mch[l]), sqr(mch[j] / mch[l]))
300 + mch[j] / mch[l] * R[a][j][l]
301 * f30(sqr(mstop[i] / mch[l]), sqr(mch[j] / mch[l])))
302 - SQRT2 * A[a][0] * T1[a][i][k] * mt_muW * mch[j] / sqr(mstop[k]) * delta(l, j)
303 * f30(sqr(mstop[i] / mstop[k]), sqr(mch[j] / mstop[k]))
304 + sqr(ma[a] / mch[j]) * delta(i, k)
305 * (U[1][j] * V[0][l]
306 * f50(sqr(mstop[i] / mch[j]), sqr(mch[l] / mch[j]), sqr(m_snutau / mch[l]))
307 - mch[l] / mch[j] * U[1][l] * V[0][j]
308 * f60(sqr(mstop[i] / mch[j]), sqr(mch[l] / mch[j]), sqr(m_snutau / mch[l])))
309 );
310 }
311 }
312 }
313 }
314 }
315 cq2_c *= -ml / 4.0 * sqr(tanb);
316
317 complex_t ca_c = 0.0;
318 for (int i = 0; i < 3; ++i) {
319 for (int j = 0; j < 2; ++j) {
320 for (int l = 0; l < 2; ++l) {
321 ca_c += I * tanb / SQRT2 * G1[i][i][j][l]
322 * (v_deltam_s * delta(l, j) * std::abs(mch[j] / mW) * f80(sqr(mstop[i] / mch[j]))
323 - (R[0][j][l] * std::abs(mch[j] / mch[l])
324 * f30(sqr(mstop[i] / mch[l]), sqr(mch[j] / mch[l]))
325 - R[0][l][j] * f40(sqr(mstop[i] / mch[l]), sqr(mch[j] / mch[l]))));
326 }
327 }
328 }
329
330 complex_t cq1 = (cq1_h + cq1_c) * mb_muW / sw2 / epsfac;
331 complex_t cq2 = (cq2_h + cq2_c) * mb_muW / sw2 / epsfac;
332 const complex_t ca = ca_h + ca_c;
333
334 // This is the first SuperIso threshold term. The lower-scale cases are
335 // handled in GroupDef_BScalar so the contribution is not counted twice.
336 if (ma[0] > Qmatch) {
337 cq2 += -v_deltam_s / 2.0 * mb_muW / sw2 * ml * ca / sqr(ma[0]);
338 }
339
340 LOG_INFO("NMSSM scalar matching enabled:",
341 "lambda=", lambda, "kappa=", kappa_nmssm,
342 "A_lambda=", A_lambda, "mu_eff=", lambda_s,
343 "M_h3=", mh[2], "M_a2=", ma[1]);
344
345 return {true, scalar_t(cq1), scalar_t(cq2)};
346}
347
348} // namespace nmssm_scalar_matching
ParameterType
ParameterType
Definition GeneralEnum.h:633
LOG_INFO
#define LOG_INFO(...)
Macro for logging informational messages.
Definition Logger.h:39
ParameterProvider.h
High-level access to parameter values and uncertainties.
complex_t
std::complex< double > complex_t
Convenience alias for std::complex<double>.
Definition Utils.h:35
ParamSrc
Lightweight view over a set of source parameters keyed by ParamId.
Definition SourcesView.h:249
ParamSrc::get_val
scalar_t get_val(const ParamId &id) const
Retrieves the current value of a parameter.
Definition SourceView.cpp:111
ParameterProvider
Provides access to parameter values, errors, and existence checks.
Definition ParameterProvider.h:53
ParameterProvider::exists
bool exists(const ParamId &pid) const
Checks if a parameter identified by ParamId exists.
Definition ParameterProvider.cpp:20
scalar_t
Definition scalar.h:11
I
constexpr std::complex< double > I
Definition constants.h:20
nmssm_scalar_matching
Definition NMSSMScalarMatching.cpp:15
nmssm_scalar_matching::compute
Result compute(const ParamSrc &src, int lepton_mass_slot)
Definition NMSSMScalarMatching.cpp:85
nmssm_scalar_matching::is_active
bool is_active()
Definition NMSSMScalarMatching.cpp:79
real
double real(const scalar_t &z)
Definition scalar.cpp:83
f40
double f40(double x, double y)
Definition special_SM.cpp:122
f30
double f30(double x, double y)
Definition special_SM.cpp:108
f60
double f60(double x, double y, double z)
Definition special_SM.cpp:159
f80
double f80(double x)
Definition special_SM.cpp:195
f50
double f50(double x, double y, double z)
Definition special_SM.cpp:136
special_SUSY.h
LhaID
Represents an identifier of a LHA element, possibly containing several sub-ids.
Definition LhaID.h:56
nmssm_scalar_matching::Result
Definition NMSSMScalarMatching.h:9