flexiblesusy is hosted by Hepforge, IPPP Durham
FlexibleSUSY
ckm.cpp
Go to the documentation of this file.
1// ====================================================================
2// This file is part of FlexibleSUSY.
3//
4// FlexibleSUSY is free software: you can redistribute it and/or modify
5// it under the terms of the GNU General Public License as published
6// by the Free Software Foundation, either version 3 of the License,
7// or (at your option) any later version.
8//
9// FlexibleSUSY is distributed in the hope that it will be useful, but
10// WITHOUT ANY WARRANTY; without even the implied warranty of
11// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12// General Public License for more details.
13//
14// You should have received a copy of the GNU General Public License
15// along with FlexibleSUSY. If not, see
16// <http://www.gnu.org/licenses/>.
17// ====================================================================
18
19#include "ckm.hpp"
20#include "error.hpp"
21#include "ew_input.hpp"
22
23#include <cmath>
24#include <limits>
25
26namespace flexiblesusy {
27
28namespace {
29
30bool is_zero(double x) noexcept
31{
32 return std::abs(x) <= std::numeric_limits<double>::epsilon();
33}
34
35double sqr(double x) { return x*x; }
36
37double pow3(double x) { return x*x*x; }
38
39double pow4(double x) { return sqr(sqr(x)); }
40
41int sign(double x) { return x >= 0.0 ? 1 : -1; }
42
43} // anonymous namespace
44
46{
47 theta_12 = 0.;
48 theta_13 = 0.;
49 theta_23 = 0.;
50 delta = 0.;
51}
52
54{
59}
60
65void CKM_parameters::set_from_wolfenstein(double lambdaW, double aCkm,
66 double rhobar, double etabar)
67{
68 if (std::abs(lambdaW) > 1.) throw SetupError("Error: Wolfenstein lambda out of range!");
69 if (std::abs(aCkm) > 1.) throw SetupError("Error: Wolfenstein A parameter out of range!");
70 if (std::abs(rhobar) > 1.) throw SetupError("Error: Wolfenstein rho-bar parameter out of range!");
71 if (std::abs(etabar) > 1.) throw SetupError("Error: Wolfenstein eta-bar parameter out of range!");
72
73 theta_12 = std::asin(lambdaW);
74 theta_23 = std::asin(aCkm * sqr(lambdaW));
75
76 const double lambdaW3 = pow3(lambdaW);
77 const double lambdaW4 = pow4(lambdaW);
78
79 const std::complex<double> rpe(rhobar, etabar);
80 const std::complex<double> V13conj = aCkm * lambdaW3 * rpe
81 * std::sqrt(1.0 - sqr(aCkm) * lambdaW4) /
82 std::sqrt(1.0 - sqr(lambdaW)) / (1.0 - sqr(aCkm) * lambdaW4 * rpe);
83
84 if (std::isfinite(std::real(V13conj)) && std::isfinite(std::imag(V13conj))) {
85 theta_13 = std::asin(std::abs(V13conj));
86 delta = std::arg(V13conj);
87 }
88}
89
94void CKM_parameters::get_wolfenstein(double& lambdaW, double& aCkm,
95 double& rhobar, double& etabar) const
96{
97 const double sin_12 = std::sin(theta_12);
98 const double sin_13 = std::sin(theta_13);
99 const double sin_23 = std::sin(theta_23);
100
101 // Eq. (11.4) from PDG
102 lambdaW = sin_12;
103 aCkm = sin_23 / sqr(lambdaW);
104
105 if (!std::isfinite(aCkm)) {
106 aCkm = 0.0;
107 }
108
109 const double c = std::sqrt((1.0 - sqr(sin_23)) / (1.0 - sqr(lambdaW)));
110 const std::complex<double> eid(std::polar(1.0, delta));
111 const std::complex<double> r(sin_13 * eid /
112 (c * aCkm * pow3(lambdaW) + sin_13 * eid * sqr(sin_23)));
113
114 rhobar = std::isfinite(std::real(r)) ? std::real(r) : 0.;
115 etabar = std::isfinite(std::imag(r)) ? std::imag(r) : 0.;
116}
117
118Eigen::Matrix<double,3,3> CKM_parameters::get_real_ckm() const
119{
120 const std::complex<double> eID(std::polar(1.0, delta));
121 const double s12 = std::sin(theta_12);
122 const double s13 = std::sin(theta_13);
123 const double s23 = std::sin(theta_23);
124 const double c12 = std::cos(theta_12);
125 const double c13 = std::cos(theta_13);
126 const double c23 = std::cos(theta_23);
127
128 // set phase factor e^(i delta) to +1 or -1 depending on the sign
129 // of its real part
130 const int pf = sign(std::real(eID));
131
132 Eigen::Matrix<double,3,3> ckm_matrix;
133 ckm_matrix(0, 0) = c12 * c13;
134 ckm_matrix(0, 1) = s12 * c13;
135 ckm_matrix(0, 2) = pf * s13;
136 ckm_matrix(1, 0) = -s12 * c23 - pf * c12 * s23 * s13;
137 ckm_matrix(1, 1) = c12 * c23 - pf * s12 * s23 * s13;
138 ckm_matrix(1, 2) = s23 * c13;
139 ckm_matrix(2, 0) = s12 * s23 - pf * c12 * c23 * s13;
140 ckm_matrix(2, 1) = -c12 * s23 - pf * s12 * c23 * s13;
141 ckm_matrix(2, 2) = c23 * c13;
142
143 return ckm_matrix;
144}
145
146Eigen::Matrix<std::complex<double>,3,3> CKM_parameters::get_complex_ckm() const
147{
148 const std::complex<double> eID(std::polar(1.0, delta));
149 const double s12 = std::sin(theta_12);
150 const double s13 = std::sin(theta_13);
151 const double s23 = std::sin(theta_23);
152 const double c12 = std::cos(theta_12);
153 const double c13 = std::cos(theta_13);
154 const double c23 = std::cos(theta_23);
155
156 Eigen::Matrix<std::complex<double>,3,3> ckm_matrix;
157 ckm_matrix(0, 0) = c12 * c13;
158 ckm_matrix(0, 1) = s12 * c13;
159 ckm_matrix(0, 2) = s13 / eID;
160 ckm_matrix(1, 0) = -s12 * c23 - c12 * s23 * s13 * eID;
161 ckm_matrix(1, 1) = c12 * c23 - s12 * s23 * s13 * eID;
162 ckm_matrix(1, 2) = s23 * c13;
163 ckm_matrix(2, 0) = s12 * s23 - c12 * c23 * s13 * eID;
164 ckm_matrix(2, 1) = -c12 * s23 - s12 * c23 * s13 * eID;
165 ckm_matrix(2, 2) = c23 * c13;
166
167 return ckm_matrix;
168}
169
170void CKM_parameters::to_pdg_convention(Eigen::Matrix<double,3,3>& Vu,
171 Eigen::Matrix<double,3,3>& Vd,
172 Eigen::Matrix<double,3,3>& Uu,
173 Eigen::Matrix<double,3,3>& Ud)
174{
175 Eigen::Matrix<double,3,3> ckm(Vu*Vd.adjoint());
176 to_pdg_convention(ckm, Vu, Vd, Uu, Ud);
177}
178
179void CKM_parameters::to_pdg_convention(Eigen::Matrix<double,3,3>& ckm,
180 Eigen::Matrix<double,3,3>& Vu,
181 Eigen::Matrix<double,3,3>& Vd,
182 Eigen::Matrix<double,3,3>& Uu,
183 Eigen::Matrix<double,3,3>& Ud)
184{
185 Eigen::Matrix<double,3,3> signs_U(Eigen::Matrix<double,3,3>::Identity()),
186 signs_D(Eigen::Matrix<double,3,3>::Identity());
187
188 // make diagonal elements positive
189 for (int i = 0; i < 3; i++) {
190 if (ckm(i, i) < 0.) {
191 signs_U(i, i) = -1.;
192 for (int j = 0; j < 3; j++) {
193 ckm(i, j) *= -1;
194 }
195 }
196 }
197
198 // make 12 element positive while keeping diagonal elements positive
199 if (ckm(0, 1) < 0.) {
200 signs_D(1, 1) = -1.;
201 signs_U(1, 1) *= -1;
202 for (int j = 0; j < 3; j++) {
203 ckm(1, j) *= -1;
204 ckm(j, 1) *= -1;
205 }
206 }
207
208 // make 23 element positive while keeping diagonal elements positive
209 if (ckm(1, 2) < 0.) {
210 signs_D(2, 2) = -1.;
211 signs_U(2, 2) *= -1;
212 for (int j = 0; j < 3; j++) {
213 ckm(2, j) *= -1;
214 ckm(j, 2) *= -1;
215 }
216 }
217
218 Vu = signs_U * Vu;
219 Uu = signs_U * Uu;
220 Vd = signs_D * Vd;
221 Ud = signs_D * Ud;
222}
223
224void CKM_parameters::to_pdg_convention(Eigen::Matrix<std::complex<double>,3,3>& Vu,
225 Eigen::Matrix<std::complex<double>,3,3>& Vd,
226 Eigen::Matrix<std::complex<double>,3,3>& Uu,
227 Eigen::Matrix<std::complex<double>,3,3>& Ud)
228{
229 Eigen::Matrix<std::complex<double>,3,3> ckm(Vu*Vd.adjoint());
230 to_pdg_convention(ckm, Vu, Vd, Uu, Ud);
231}
232
233namespace {
234
235template<typename T>
236std::complex<T> phase(const std::complex<T>& z) noexcept
237{
238 T r = std::abs(z);
239 return is_zero(r) ? 1 : z/r;
240}
241
243(const Eigen::Matrix<std::complex<double>,3,3>& ckm,
244 const std::complex<double>& p,
245 std::complex<double>& o,
246 Eigen::DiagonalMatrix<std::complex<double>,3>& l,
247 Eigen::DiagonalMatrix<std::complex<double>,3>& r)
248{
249 o = std::conj(phase(p * ckm(0,2)));
250 l.diagonal().bottomRightCorner<2,1>() = (o * ckm.bottomRightCorner<2,1>()).
251 unaryExpr([] (const std::complex<double>& z) { return phase<double>(z); }).conjugate();
252 r.diagonal().topLeftCorner<2,1>() = (o * ckm.topLeftCorner<1,2>()).
253 unaryExpr([] (const std::complex<double>& z) { return phase<double>(z); }).adjoint();
254}
255
257double sanitize_hypot(double sc) noexcept
258{
259 if (sc < -1.) { sc = -1.0; }
260 if (sc > 1.) { sc = 1.0; }
261 return sc;
262}
263
264} // anonymous namespace
265
266void CKM_parameters::to_pdg_convention(Eigen::Matrix<std::complex<double>,3,3>& ckm,
267 Eigen::Matrix<std::complex<double>,3,3>& Vu,
268 Eigen::Matrix<std::complex<double>,3,3>& Vd,
269 Eigen::Matrix<std::complex<double>,3,3>& Uu,
270 Eigen::Matrix<std::complex<double>,3,3>& Ud)
271{
272 std::complex<double> o;
273 Eigen::DiagonalMatrix<std::complex<double>,3> l(1,1,1), r(1,1,1);
274
275 const double s13 = sanitize_hypot(std::abs(ckm(0,2)));
276 const double c13_sq = 1. - sqr(s13);
277
278 if (is_zero(c13_sq)) {
279 o = std::conj(phase(ckm(0,2)));
280 r.diagonal().block<2,1>(0,0) = (o * ckm.block<1,2>(1,0)).unaryExpr(
281 [] (const std::complex<double>& z) { return phase<double>(z); }).adjoint();
282 l.diagonal()[2] = std::conj(phase(o * r.diagonal()[1] * ckm(2,1)));
283 } else {
284 const double c13 = std::sqrt(c13_sq);
285 const double s12 = sanitize_hypot(std::abs(ckm(0,1)) / c13);
286 const double c12 = std::sqrt(1 - sqr(s12));
287 const double s23 = sanitize_hypot(std::abs(ckm(1,2)) / c13);
288 const double c23 = std::sqrt(1 - sqr(s23));
289 const double side1 = s12*s23;
290 const double side2 = c12*c23*s13;
291 const double cosdelta = sanitize_hypot(
292 is_zero(side1*side2) ? 1 // delta is removable
293 : (sqr(side1)+sqr(side2)-std::norm(ckm(2,0))) / (2*side1*side2));
294 const double sindelta = std::sqrt(1 - sqr(cosdelta));
295 const std::complex<double> p(cosdelta, sindelta); // Exp[I delta]
296 calc_phase_factors(ckm, p, o, l, r);
297 const Eigen::Array<double,2,2>
298 imagBL{(o * l * ckm * r).bottomLeftCorner<2,2>().imag()};
299 if (!((imagBL <= 0).all() || (imagBL >= 0).all())) {
300 calc_phase_factors(ckm, std::conj(p), o, l, r);
301 }
302 }
303
304 Vu.transpose() *= l * o;
305 Vd.transpose() *= r.diagonal().conjugate().asDiagonal();
306 Uu.transpose() *= (l * o).diagonal().conjugate().asDiagonal();
307 Ud.transpose() *= r;
308 ckm = Vu * Vd.adjoint();
309}
310
311} // namespace flexiblesusy
Spectrum generator was not setup correctly.
Definition: error.hpp:46
constexpr double CKM_THETA23
From Vcb/Vtb in global CKM fit, PDG.
Definition: ew_input.hpp:78
constexpr double CKM_THETA13
From Vub in global CKM fit, PDG.
Definition: ew_input.hpp:77
constexpr double CKM_THETA12
From Vus/Vud in global CKM fit, PDG.
Definition: ew_input.hpp:76
double sanitize_hypot(double sc) noexcept
restrict sin or cos to interval [-1,1]
Definition: ckm.cpp:257
void calc_phase_factors(const Eigen::Matrix< std::complex< double >, 3, 3 > &ckm, const std::complex< double > &p, std::complex< double > &o, Eigen::DiagonalMatrix< std::complex< double >, 3 > &l, Eigen::DiagonalMatrix< std::complex< double >, 3 > &r)
Definition: ckm.cpp:243
std::complex< T > phase(const std::complex< T > &z) noexcept
Definition: ckm.cpp:236
constexpr T norm(const Complex< T > &z) noexcept
Definition: complex.hpp:66
T sqr(T x)
Definition: mathdefs.hpp:62
const Real epsilon
Definition: mathdefs.hpp:18
constexpr T arg(const Complex< T > &z) noexcept
Definition: complex.hpp:42
std::enable_if_t< std::is_unsigned< T >::value, bool > is_zero(T a, T prec=std::numeric_limits< T >::epsilon()) noexcept
compares a number for being close to zero
Definition: numerics2.hpp:46
T sign(T x)
Definition: mathdefs.hpp:65
constexpr Complex< T > conj(const Complex< T > &z) noexcept
Definition: complex.hpp:48
void set_from_wolfenstein(double, double, double, double)
Definition: ckm.cpp:65
Eigen::Matrix< double, 3, 3 > get_real_ckm() const
Definition: ckm.cpp:118
static void to_pdg_convention(Eigen::Matrix< double, 3, 3 > &, Eigen::Matrix< double, 3, 3 > &, Eigen::Matrix< double, 3, 3 > &, Eigen::Matrix< double, 3, 3 > &, Eigen::Matrix< double, 3, 3 > &)
Definition: ckm.cpp:179
Eigen::Matrix< std::complex< double >, 3, 3 > get_complex_ckm() const
Definition: ckm.cpp:146
void get_wolfenstein(double &, double &, double &, double &) const
Definition: ckm.cpp:94