47 auto zero() {
return std::complex<T>(0.0,0.0); }
49 auto one() {
return std::complex<T>(1.0,0.0); }
52 static const double kK1 = 2.53386551601e-00;
53 static const double kK2 = 4.62711492217e-09;
61 this->SetMix(0.,0.,0.,0.);
62 this->SetDeltaMsqrs(0.,0.);
70 const T& dms21,
const T& dms32)
72 this->SetMix(th12, th23, th13, deltacp);
73 this->SetDeltaMsqrs(dms21, dms32);
83 <<
"| "<<M[0][0]<<
"\t"<<M[0][1]<<
"\t"<<M[0][2]<<
" |\n" 84 <<
"| "<<M[1][0]<<
"\t"<<M[1][1]<<
"\t"<<M[1][2]<<
" |\n" 85 <<
"| "<<M[2][0]<<
"\t"<<M[2][1]<<
"\t"<<M[2][2]<<
" |\n" 100 T s12, s23, s13, c12, c23, c13;
103 s12 =
sin(th12); s23 =
sin(th23); s13 =
sin(th13);
104 c12 =
cos(th12); c23 =
cos(th23); c13 =
cos(th13);
108 fU[0][0] =
complex(c12*c13, 0);
109 fU[0][1] =
complex(s12*c13, 0);
110 fU[0][2] = s13*
exp(-idelta);
112 fU[1][0] = -s12*c23-c12*s23*s13*
exp(idelta);
113 fU[1][1] = c12*c23-s12*s23*s13*
exp(idelta);
114 fU[1][2] =
complex(s23*c13, 0);
116 fU[2][0] = s12*s23-c12*c23*s13*
exp(idelta);
117 fU[2][1] = -c12*s23-s12*c23*s13*
exp(idelta);
118 fU[2][2] =
complex(c23*c13, 0);
121 for (i=0; i<3; ++
i) {
122 for (j=0; j<3; ++
j) {
123 fUtran[
i][
j] = fU[
j][
i];
124 fUstar[
i][
j] = conj(fU[i][j]);
125 fUdagg[
i][
j] = conj(fU[j][i]);
138 template <
typename T>
142 double s1, s2, s3,
c1,
c2,
c3;
143 std::complex<double>
id(0.0,d);
144 s1 =
sin(th1); s2 =
sin(th2); s3 =
sin(th3);
145 c1 =
cos(th1); c2 =
cos(th2); c3 =
cos(th3);
152 fU[1][1] = c1*c2*c3+s2*s3*
exp(
id);
153 fU[1][2] = c1*c2*s3-s2*c3*
exp(
id);
156 fU[2][1] = c1*s2*c3-c2*s3*
exp(
id);
157 fU[2][2] = c1*s2*s3+c2*c3*
exp(
id);
160 for (i=0; i<3; ++
i) {
161 for (j=0; j<3; ++
j) {
162 fUtran[
i][
j] = fU[
j][
i];
163 fUstar[
i][
j] = conj(fU[i][j]);
164 fUdagg[
i][
j] = conj(fU[j][i]);
171 template <
typename T>
175 <<
"|"<<fdmsqr[0][0]<<
"\t"<<fdmsqr[0][1]<<
"\t"<<fdmsqr[0][2]<<
"|\n" 176 <<
"|"<<fdmsqr[1][0]<<
"\t"<<fdmsqr[1][1]<<
"\t"<<fdmsqr[1][2]<<
"|\n" 177 <<
"|"<<fdmsqr[2][0]<<
"\t"<<fdmsqr[2][1]<<
"\t"<<fdmsqr[2][2]<<
"|" 186 template <
typename T>
198 if (dm21==0.0) {msqr[0] -= 0.5*eps; msqr[1] += 0.5*eps; }
199 if (dm32==0.0) {msqr[1] -= 0.5*eps; msqr[2] += 0.5*eps; }
203 for (
int i=0;
i<3; ++
i) {
204 for (
int j=0;
j<3; ++
j) {
210 fdmsqr[
i][
j] = -(msqr[
i] - msqr[
j]);
219 template <
typename T>
223 for (i=0; i<3; ++
i) {
224 for (j=0; j<3; ++
j) {
226 for (k=0; k<3; ++k) {
227 A[
i][
j] += B[
i][k]*C[k][
j];
240 template <
typename T>
249 for (a=0; a<3; ++
a) {
250 for (b=0; b<3; ++
b) {
252 for (i=0; i<3; ++
i) {
254 A[
a][
b] += U[
a][
i]*
exp(phase)*Udagg[
i][
b];
264 template <
typename T>
275 for (k=0; k<3; ++k) {
276 for (j=0; j<3; ++
j) {
277 twoEH[k][
j] = zero<T>();
278 if (k==j) twoEH[k][
j] =
complex(dmsqr[j][0], 0);
279 twoEH[k][
j] -= k2r2GNeE*U[0][
j]*Udagg[k][0];
291 template <
typename T>
298 this->Multi(tmp, X, Udagg);
299 this->Multi(A, U, tmp);
307 template <
typename T>
316 static const double One[3][3] = {{1.,0.,0.},
331 for (a=0; a<3; ++
a) {
332 for (b=0; b<3; ++
b) {
333 EHM0[
a][
b] = twoEH[
a][
b]-Msqr[0]*One[
a][
b];
334 EHM1[
a][
b] = twoEH[
a][
b]-Msqr[1]*One[
a][
b];
335 EHM2[
a][
b] = twoEH[
a][
b]-Msqr[2]*One[
a][
b];
338 this->Multi(EHM21,EHM2,EHM1);
339 this->Multi(EHM20,EHM2,EHM0);
340 this->Multi(EHM10,EHM1,EHM0);
343 for (a=0; a<3; ++
a) {
344 for (b=0; b<3; ++
b) {
346 for (k=0; k<3; ++k) {
349 X[
a][
b] += (EHM21[
a][
b]/(dMsqr[k][2]*dMsqr[k][1]))*phase;
352 X[
a][
b] += (EHM20[
a][
b]/(dMsqr[k][2]*dMsqr[k][0]))*phase;
355 X[
a][
b] += (EHM10[
a][
b]/(dMsqr[k][1]*dMsqr[k][0]))*phase;
367 template <
typename T>
373 static const auto k2PiOver3 = 2.0*
M_PI/3.0;
374 auto alpha2 = alpha*alpha;
375 auto alpha3 = alpha*alpha2;
376 auto alpha2Minus3beta = alpha2-3.0*
beta;
380 (2.0*alpha3 - 9.0*alpha*beta + 27.0*gamma)/
381 (2.0*
pow(alpha2Minus3beta,1.5));
387 if (
fabs(arg-1.0)>1.E-9) abort();
388 if (arg<-1.0) arg = -1.00;
389 if (arg>+1.0) arg = +1.00;
394 auto theta0 =
acos(arg)/3.0;
395 auto theta1 = theta0-k2PiOver3;
396 auto theta2 = theta0+k2PiOver3;
399 auto fac = -2.0*
sqrt(alpha2Minus3beta)/3.0;
400 auto constant = -alpha/3.0;
403 Msqr[0] = fac*
cos(theta0) + constant;
404 Msqr[1] = fac*
cos(theta1) + constant;
405 Msqr[2] = fac*
cos(theta2) + constant;
413 template <
typename T>
426 alpha = k2r2EGNe + dmsqr[0][1] + dmsqr[0][2];
429 dmsqr[0][1]*dmsqr[0][2] +
430 k2r2EGNe*(dmsqr[0][1]*(1.0-
norm(U[0][1])) +
431 dmsqr[0][2]*(1.0-
norm(U[0][2])));
433 gamma = k2r2EGNe*dmsqr[0][1]*dmsqr[0][2]*
norm(U[0][0]);
440 template <
typename T>
448 T MsqrTmp[3] = {-99.9E9,-99.9E9,-99.9E9};
449 int flg[3] = {0,0,0};
453 for (i=0; i<3; ++
i) {
456 for (j=0; j<3; ++
j) {
458 auto delta =
abs(MsqrVac[i] - dmsqr[j][0]);
461 if (best>1.
E-9) abort();
462 if (k<0 || k>2) abort();
464 MsqrTmp[
i] = Msqr[k];
467 for (i=0; i<3; ++
i)
if (flg[i]!=1) abort();
469 for (i=0; i<3; ++
i) {
470 Msqr[
i] = MsqrTmp[
i];
471 for (j=0; j<3; ++
j) {
472 dMsqr[
i][
j] = (MsqrTmp[
i] - MsqrTmp[
j]);
481 template <
typename T>
488 if (anti>0) this->EvalEqn2(A, fU, fUdagg, fdmsqr, L, E);
489 else if (anti<0) this->EvalEqn2(A, fUstar, fUtran, fdmsqr, L, E);
491 this->Multi(Aout, A, fA);
492 for (i=0; i<3; ++
i) {
493 for (j=0; j<3; ++
j) {
494 fA[
i][
j] = Aout[
i][
j];
503 template <
typename T>
506 static const double Gf = 1.166371E-5;
513 T alpha,
beta, gamma;
514 T alpha0, beta0, gamma0;
521 this->EvalEqn5(twoEH, fU, fUdagg, fdmsqr, E, Gf, Ne);
524 this->EvalEqn22(alpha, beta, gamma, E, Gf, Ne, fdmsqr, fU);
525 this->EvalEqn21(Msqr, alpha, beta, gamma);
529 this->EvalEqn22(alpha0, beta0, gamma0, E, 0.0, 0.0, fdmsqr, fU);
530 this->EvalEqn21(MsqrV, alpha0, beta0, gamma0);
531 this->SortEigenvalues(dMsqr, fdmsqr, MsqrV, Msqr);
534 this->EvalEqn11(X, L, E, twoEH, Msqr, dMsqr);
535 this->EvalEqn10(A, fU, X, fUdagg);
540 this->EvalEqn5(twoEH, fUstar, fUtran, fdmsqr, E, -Gf, Ne);
541 this->EvalEqn22(alpha, beta, gamma, E, -Gf, Ne, fdmsqr, fUstar);
542 this->EvalEqn21(Msqr, alpha, beta, gamma);
543 this->EvalEqn22(alpha0, beta0, gamma0, E, 0.0, 0.0, fdmsqr, fUstar);
544 this->EvalEqn21(MsqrV, alpha0, beta0, gamma0);
545 this->SortEigenvalues(dMsqr, fdmsqr, MsqrV, Msqr);
546 this->EvalEqn11(X, L, E, twoEH, Msqr, dMsqr);
547 this->EvalEqn10(A, fUstar, X, fUtran);
552 this->Multi(Aout, A, fA);
553 for (i=0; i<3; ++
i) {
554 for (j=0; j<3; ++
j) {
555 fA[
i][
j] = Aout[
i][
j];
564 template <
typename T>
567 const std::list<double>& Ne,
570 if (L.size()!=Ne.size()) abort();
573 auto Ni (Ne.begin());
574 for (; Li!=Lend; ++Li, ++Ni) {
576 static const double kRhoCutoff = 1.0E-6;
577 if (*Ni<kRhoCutoff) this->PropVacuum(*Li, E, anti);
578 else this->PropMatter(*Li, E, *Ni, anti);
587 template <
typename T>
591 for (i=0; i<3; ++
i) {
592 for (j=0; j<3; ++
j) {
593 if (i==j) fA[
i][
j] = one<T>();
594 else fA[
i][
j] = zero<T>();
605 template <
typename T>
610 return norm(fA[i][j]);
::xsd::cxx::tree::id< char, ncname > id
void EvalEqn2(complex A[][3], const complex U[][3], const complex Udagg[][3], const T dmsqr[][3], double L, double E)
fvar< T > fabs(const fvar< T > &x)
static const double kK2
mole/GeV^2/cm^3 to eV
void EvalEqn11(complex X[][3], double L, double E, const complex twoEH[][3], const T Msqr[], const T dMsqr[][3])
void PrintDeltaMsqrs() const
Print the matrix of mass-squared differneces.
void SortEigenvalues(T dMsqr[][3], const T dmsqr[][3], const T MsqrVac[], T Msqr[])
static const double kGeV2eV
GeV to eV.
static constexpr double L
void EvalEqn10(complex A[][3], const complex U[][3], const complex X[][3], const complex Udagg[][3])
void EvalEqn22(T &alpha, T &beta, T &gamma, double E, double Gf, double Ne, const T dmsqr[][3], const complex U[][3])
fvar< T > exp(const fvar< T > &x)
void PropMatter(double L, double E, double Ne, int anti)
void PropVacuum(double L, double E, int anti)
void SetMixBWCP(double th1, double th2, double th3, double deltacp)
Oscillation probability calculators.
void Multi(complex A[][3], const complex B[][3], const complex C[][3])
assert(nhit_max >=nhit_nbins)
void SetMix(const T &th12, const T &th23, const T &th13, const T &deltacp)
void SetDeltaMsqrs(const T &dm21, const T &dm32)
void PrintMix() const
Print the oscillation matrix.
void EvalEqn5(complex twoEH[][3], const complex U[][3], const complex Udagg[][3], const T dmsqr[][3], double E, double Gf, double Ne)
static const double kK1
(1/2)*(1000/hbarc)
void EvalEqn21(T Msqr[], const T &alpha, const T &beta, const T &gamma)
void DumpMatrix(const complex M[][3]) const
Utility to print a genetrix complex matrix M.