35double sqr(
double x) {
return x*x; }
37double pow3(
double x) {
return x*x*x; }
41int sign(
double x) {
return x >= 0.0 ? 1 : -1; }
66 double rhobar,
double etabar)
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!");
76 const double lambdaW3 =
pow3(lambdaW);
77 const double lambdaW4 =
pow4(lambdaW);
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);
84 if (std::isfinite(std::real(V13conj)) && std::isfinite(std::imag(V13conj))) {
85 theta_13 = std::asin(std::abs(V13conj));
95 double& rhobar,
double& etabar)
const
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);
103 aCkm = sin_23 /
sqr(lambdaW);
105 if (!std::isfinite(aCkm)) {
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)));
114 rhobar = std::isfinite(std::real(r)) ? std::real(r) : 0.;
115 etabar = std::isfinite(std::imag(r)) ? std::imag(r) : 0.;
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);
130 const int pf =
sign(std::real(eID));
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;
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);
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;
171 Eigen::Matrix<double,3,3>& Vd,
172 Eigen::Matrix<double,3,3>& Uu,
173 Eigen::Matrix<double,3,3>& Ud)
175 Eigen::Matrix<double,3,3> ckm(Vu*Vd.adjoint());
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)
185 Eigen::Matrix<double,3,3> signs_U(Eigen::Matrix<double,3,3>::Identity()),
186 signs_D(Eigen::Matrix<double,3,3>::Identity());
189 for (
int i = 0; i < 3; i++) {
190 if (ckm(i, i) < 0.) {
192 for (
int j = 0; j < 3; j++) {
199 if (ckm(0, 1) < 0.) {
202 for (
int j = 0; j < 3; j++) {
209 if (ckm(1, 2) < 0.) {
212 for (
int j = 0; j < 3; j++) {
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)
229 Eigen::Matrix<std::complex<double>,3,3> ckm(Vu*Vd.adjoint());
236std::complex<T>
phase(
const std::complex<T>& z)
noexcept
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)
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();
259 if (sc < -1.) { sc = -1.0; }
260 if (sc > 1.) { sc = 1.0; }
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)
272 std::complex<double> o;
273 Eigen::DiagonalMatrix<std::complex<double>,3> l(1,1,1), r(1,1,1);
276 const double c13_sq = 1. -
sqr(s13);
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)));
284 const double c13 = std::sqrt(c13_sq);
286 const double c12 = std::sqrt(1 -
sqr(s12));
288 const double c23 = std::sqrt(1 -
sqr(s23));
289 const double side1 = s12*s23;
290 const double side2 = c12*c23*s13;
294 const double sindelta = std::sqrt(1 -
sqr(cosdelta));
295 const std::complex<double> p(cosdelta, sindelta);
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())) {
304 Vu.transpose() *= l * o;
305 Vd.transpose() *= r.diagonal().conjugate().asDiagonal();
306 Uu.transpose() *= (l * o).diagonal().conjugate().asDiagonal();
308 ckm = Vu * Vd.adjoint();
Spectrum generator was not setup correctly.
double sanitize_hypot(double sc) noexcept
restrict sin or cos to interval [-1,1]
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)
std::complex< T > phase(const std::complex< T > &z) noexcept
constexpr T norm(const Complex< T > &z) noexcept
constexpr T arg(const Complex< T > &z) noexcept
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
constexpr Complex< T > conj(const Complex< T > &z) noexcept
void set_from_wolfenstein(double, double, double, double)
void reset_to_observation()
Eigen::Matrix< double, 3, 3 > get_real_ckm() const
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 > &)
Eigen::Matrix< std::complex< double >, 3, 3 > get_complex_ckm() const
void get_wolfenstein(double &, double &, double &, double &) const