PMNS.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // $Id: PMNS.cxx,v 1.4 2012-09-20 21:42:52 greenc Exp $
3 //
4 // Implementation of oscillations of neutrinos in matter in an
5 // three-neutrino framework based on the following reference:
6 //
7 //......................................................................
8 //
9 // PHYSICAL REVIEW D VOLUME 22, NUMBER 11 1 DECEMBER 1980
10 //
11 // Matter effects on three-neutrino oscillation
12 //
13 // V. Barger and K. Whisnant
14 // Physics Department, U. of Wisconsin, Madison, Wisconsin 53706
15 //
16 // S. Pakvasa
17 // Physics Department, U. of Hawaii at Manoa, Honolulu, Hawaii 96822
18 //
19 // R.J.N. Phillips
20 // Rutherford Laboratory, Chilton, Didcot, Oxon, England
21 // (Received 4 August 1980)
22 //
23 // 22 2718
24 // --
25 //......................................................................
26 //
27 // Throughout I have taken:
28 // - L to be the neutrino flight distance in km
29 // - E to be the neutrino energy in GeV
30 // - dmsqr to be the differences between the mass-squares in eV^2
31 // - Ne to be the electron number density in mole/cm^3
32 // - theta12,theta23,theta13,deltaCP to be in radians
33 //
34 // messier@indiana.edu
35 ////////////////////////////////////////////////////////////////////////
36 #include "OscLib/PMNS.h"
37 
38 #include <cstdlib>
39 #include <iostream>
40 #include <cassert>
41 
42 using namespace osc;
43 
44 
45 // Some useful complex numbers
46 template<typename T>
47 auto zero() { return std::complex<T>(0.0,0.0); }
48 template<typename T>
49 auto one() { return std::complex<T>(1.0,0.0); }
50 
51 // Unit conversion constants
52 static const double kK1 = 2.53386551601e-00; ///< (1/2)*(1000/hbarc)
53 static const double kK2 = 4.62711492217e-09; ///< mole/GeV^2/cm^3 to eV
54 static const double kGeV2eV = 1.0E9; ///< GeV to eV
55 
56 //......................................................................
57 
58 template <typename T>
60 {
61  this->SetMix(0.,0.,0.,0.);
62  this->SetDeltaMsqrs(0.,0.);
63  this->Reset();
64 }
65 
66 //......................................................................
67 
68 template <typename T>
69 _PMNS<T>::_PMNS(const T& th12, const T& th23, const T& th13, const T& deltacp,
70  const T& dms21, const T& dms32)
71 {
72  this->SetMix(th12, th23, th13, deltacp);
73  this->SetDeltaMsqrs(dms21, dms32);
74  this->Reset();
75 }
76 
77 //......................................................................
78 
79 template <typename T>
80 void _PMNS<T>::DumpMatrix(const complex M[][3]) const
81 {
82  std::cout
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"
86  <<std::endl;
87 }
88 
89 //......................................................................
90 
91 template <typename T>
92 void _PMNS<T>::PrintMix() const { this->DumpMatrix(fU); }
93 
94 //......................................................................
95 
96 template <typename T>
97 void _PMNS<T>::SetMix(const T& th12, const T& th23, const T& th13, const T& deltacp)
98 {
99  int i, j;
100  T s12, s23, s13, c12, c23, c13;
101  complex idelta(0.0,deltacp);
102 
103  s12 = sin(th12); s23 = sin(th23); s13 = sin(th13);
104  c12 = cos(th12); c23 = cos(th23); c13 = cos(th13);
105 
106  // explicitly initializing complex part
107  // to circumvent problem described at top of this file
108  fU[0][0] = complex(c12*c13, 0);
109  fU[0][1] = complex(s12*c13, 0);
110  fU[0][2] = s13*exp(-idelta);
111 
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);
115 
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);
119 
120  // Compute derived forms of the mixing matrix
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]);
126  }
127  }
128 }
129 
130 ///.....................................................................
131 ///
132 /// Initialize the mixing matrix using the older form referenced by
133 /// the Barger et al paper
134 ///
135 /// \warning This should not be used except for testing. Use SetMix
136 /// above.
137 ///
138 template <typename T>
139 void _PMNS<T>::SetMixBWCP(double th1, double th2, double th3, double d)
140 {
141  int i, j;
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);
146 
147  fU[0][0] = c1;
148  fU[0][1] = s1*c3;
149  fU[0][2] = s1*s3;
150 
151  fU[1][0] = -s1*c2;
152  fU[1][1] = c1*c2*c3+s2*s3*exp(id);
153  fU[1][2] = c1*c2*s3-s2*c3*exp(id);
154 
155  fU[2][0] = -s1*s2;
156  fU[2][1] = c1*s2*c3-c2*s3*exp(id);
157  fU[2][2] = c1*s2*s3+c2*c3*exp(id);
158 
159  // Compute derived forms of the mixing matrix
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]);
165  }
166  }
167 }
168 
169 //......................................................................
170 
171 template <typename T>
173 {
174  std::cout
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]<<"|"
178  << std::endl;
179 }
180 
181 //......................................................................
182 ///
183 /// Set the mass-splittings. These are m_2^2-m_1^2 and m_3^2-m_2^2 in
184 /// eV^2
185 ///
186 template <typename T>
187 void _PMNS<T>::SetDeltaMsqrs(const T& dm21, const T& dm32)
188 {
189  double eps = 5.0E-9;
190  T msqr[3];
191 
192  msqr[0] = 0.0;
193  msqr[1] = dm21;
194  msqr[2] = dm21+dm32;
195 
196  // Degeneracies cause problems with diagonalization, so break them
197  // ever so slightly
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; }
200 
201  // Assign the mass splittings fdmsqr_ij = msqr_i - msqr_j by
202  // convention
203  for (int i=0; i<3; ++i) {
204  for (int j=0; j<3; ++j) {
205  // A somewhat subtle point: Barger et al. refer to the sign of
206  // m1-m2 being positive which implies dm^2_12>0 and
207  // dm^2_21<0. The labeling in more common use is to assume m3 is
208  // heaviest such that dm_12<0 and dm_21>0. Rather than reverse
209  // all the indices in all the equations, flip the signs here.
210  fdmsqr[i][j] = -(msqr[i] - msqr[j]);
211  }
212  }
213 }
214 
215 //......................................................................
216 ///
217 /// Compute matrix multiplication A = B*C
218 ///
219 template <typename T>
220 void _PMNS<T>::Multi(complex A[][3], const complex B[][3], const complex C[][3])
221 {
222  int i, j, k;
223  for (i=0; i<3; ++i) {
224  for (j=0; j<3; ++j) {
225  A[i][j] = zero<T>();
226  for (k=0; k<3; ++k) {
227  A[i][j] += B[i][k]*C[k][j];
228  }
229  }
230  }
231 }
232 
233 //......................................................................
234 ///
235 /// Compute Equation 2, transition matrix for propagation through
236 /// vaccum, from Barger et al:
237 ///
238 /// A(nua->nub) = sum_i U_ai exp(-1/2 i m_i^2 L/E) Udagg_ib
239 ///
240 template <typename T>
242  const complex U[][3],
243  const complex Udagg[][3],
244  const T dmsqr[][3],
245  double L,
246  double E)
247 {
248  int a, b, i;
249  for (a=0; a<3; ++a) {
250  for (b=0; b<3; ++b) {
251  A[a][b] = zero<T>();
252  for (i=0; i<3; ++i) {
253  complex phase(0.0,-kK1*dmsqr[i][0]*L/E);
254  A[a][b] += U[a][i]*exp(phase)*Udagg[i][b];
255  }
256  }
257  }
258 }
259 
260 //......................................................................
261 ///
262 /// Compute Eqn. 5, the matter 2*E*Hamiltonian in the mass basis
263 ///
264 template <typename T>
265 void _PMNS<T>::EvalEqn5(complex twoEH[][3],
266  const complex U[][3],
267  const complex Udagg[][3],
268  const T dmsqr[][3],
269  double E,
270  double Gf,
271  double Ne)
272 {
273  int j, k;
274  T k2r2GNeE = kK2*2.0*M_SQRT2*Gf*Ne*(kGeV2eV*E);
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];
280  }
281  }
282 }
283 
284 //......................................................................
285 ///
286 /// Compute Equation 10, the transition matrix for propagation across
287 /// matter slab of constant density, from Barger et al:
288 ///
289 /// A = U * X * U^(dagger) for neutrinos
290 ///
291 template <typename T>
293  const complex U[][3],
294  const complex X[][3],
295  const complex Udagg[][3])
296 {
297  complex tmp[3][3];
298  this->Multi(tmp, X, Udagg);
299  this->Multi(A, U, tmp);
300 }
301 
302 //......................................................................
303 ///
304 /// Evaluate Eqn. 11, the Lagragne formula for the matrix e(-iHt),
305 /// from Barger et al.
306 ///
307 template <typename T>
309  double L,
310  double E,
311  const complex twoEH[][3],
312  const T Msqr[],
313  const T dMsqr[][3])
314 {
315  // The identity matrix
316  static const double One[3][3] = {{1.,0.,0.},
317  {0.,1.,0.},
318  {0.,0.,1.}
319  };
320  int a, b, k;
321  complex phase;
322  complex EHM0[3][3];
323  complex EHM1[3][3];
324  complex EHM2[3][3];
325  complex EHM21[3][3];
326  complex EHM20[3][3];
327  complex EHM10[3][3];
328 
329  // There are three matrices which can apper inside the product on
330  // j!=k. Calculate them before hand
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];
336  }
337  }
338  this->Multi(EHM21,EHM2,EHM1);
339  this->Multi(EHM20,EHM2,EHM0);
340  this->Multi(EHM10,EHM1,EHM0);
341 
342  // Refer to Msqr_j as dMsqr[j][0] since only mass differences matter
343  for (a=0; a<3; ++a) {
344  for (b=0; b<3; ++b) {
345  X[a][b] = zero<T>();
346  for (k=0; k<3; ++k) {
347  phase = exp(complex(0.0,-kK1*dMsqr[k][0]*L/E));
348  if (k==0) {
349  X[a][b] += (EHM21[a][b]/(dMsqr[k][2]*dMsqr[k][1]))*phase;
350  }
351  else if (k==1) {
352  X[a][b] += (EHM20[a][b]/(dMsqr[k][2]*dMsqr[k][0]))*phase;
353  }
354  else if (k==2) {
355  X[a][b] += (EHM10[a][b]/(dMsqr[k][1]*dMsqr[k][0]))*phase;
356  } // Switch for product on j!=k
357  } // Sum on k
358  } // Loop on b
359  } // Loop on a
360 }
361 
362 //......................................................................
363 //
364 // Find the matter eigenvalues Msqr given the variables found in
365 // Eqn.22. This is Eqn.21 from Barger et. al.
366 //
367 template <typename T>
368 void _PMNS<T>::EvalEqn21(T Msqr[],
369  const T& alpha,
370  const T& beta,
371  const T& gamma)
372 {
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;
377 
378  // Argument of the acos()
379  auto arg =
380  (2.0*alpha3 - 9.0*alpha*beta + 27.0*gamma)/
381  (2.0*pow(alpha2Minus3beta,1.5));
382 
383  // Occasionally round off errors mean that arg wanders outside of
384  // its allowed range. If its way off (1 part per billion), stop the
385  // program. Otherwise, set it to its real value.
386  if (fabs(arg)>1.0) {
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;
390  }
391 
392  // The three roots, find the first by computing the acos() the
393  // others are nearby
394  auto theta0 = acos(arg)/3.0;
395  auto theta1 = theta0-k2PiOver3;
396  auto theta2 = theta0+k2PiOver3;
397 
398  // The multiplier and offset
399  auto fac = -2.0*sqrt(alpha2Minus3beta)/3.0; // Factor in front of cos() terms
400  auto constant = -alpha/3.0; // The constant offset m1^2 is irrelevant
401 
402  // And the eigenvalues themselves
403  Msqr[0] = fac*cos(theta0) + constant;
404  Msqr[1] = fac*cos(theta1) + constant;
405  Msqr[2] = fac*cos(theta2) + constant;
406 }
407 
408 ///.....................................................................
409 ///
410 /// Compute the values of the simplifying expressions alpha, beta, and
411 /// gamma. This is Eqn22 from the Barger et al paper
412 ///
413 template <typename T>
414 void _PMNS<T>::EvalEqn22(T& alpha,
415  T& beta,
416  T& gamma,
417  double E,
418  double Gf,
419  double Ne,
420  const T dmsqr[][3],
421  const complex U[][3])
422 {
423  // 2*sqrt(2)*Gf*Ne*E in units of eV^2
424  auto k2r2EGNe = kK2*2.0*M_SQRT2*Gf*Ne*(kGeV2eV*E);
425 
426  alpha = k2r2EGNe + dmsqr[0][1] + dmsqr[0][2];
427 
428  beta =
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])));
432 
433  gamma = k2r2EGNe*dmsqr[0][1]*dmsqr[0][2]*norm(U[0][0]);
434 }
435 
436 //......................................................................
437 ///
438 /// Sort out the eigenvalues
439 ///
440 template <typename T>
441 void _PMNS<T>::SortEigenvalues(T dMsqr[][3],
442  const T dmsqr[][3],
443  const T MsqrVac[],
444  T Msqr[])
445 {
446  int i, j, k;
447  T best;
448  T MsqrTmp[3] = {-99.9E9,-99.9E9,-99.9E9};
449  int flg[3] = {0,0,0};
450 
451  // Attempt to figure out which of the eigenvalues match between
452  // dmsqr and MsqrVac
453  for (i=0; i<3; ++i) {
454  best = 1.E30;
455  k = -1;
456  for (j=0; j<3; ++j) {
457  using namespace std;
458  auto delta = abs(MsqrVac[i] - dmsqr[j][0]);
459  if (delta<best) { best = delta; k = j; }
460  }
461  if (best>1.E-9) abort();
462  if (k<0 || k>2) abort();
463  flg[k] = 1;
464  MsqrTmp[i] = Msqr[k];
465  }
466  // Check that each eigenvalue is used
467  for (i=0; i<3; ++i) if (flg[i]!=1) abort();
468 
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]);
473  }
474  }
475 }
476 
477 ///.....................................................................
478 ///
479 /// Update the transition matrix for a step across the vacuum
480 ///
481 template <typename T>
482 void _PMNS<T>::PropVacuum(double L, double E, int anti)
483 {
484  int i, j;
485  complex A[3][3];
486  complex Aout[3][3];
487 
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);
490  else abort();
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];
495  }
496  }
497 }
498 
499 ///.....................................................................
500 ///
501 /// Update the transition matrix for a step across a slab of matter
502 ///
503 template <typename T>
504 void _PMNS<T>::PropMatter(double L, double E, double Ne, int anti)
505 {
506  static const double Gf = 1.166371E-5; // G_F in units of GeV^-2
507  int i, j;
508  complex twoEH[3][3];
509  complex X[3][3];
510  T Msqr[3];
511  T MsqrV[3];
512  T dMsqr[3][3];
513  T alpha, beta, gamma;
514  T alpha0, beta0, gamma0;
515  complex A[3][3];
516  complex Aout[3][3];
517 
518  // Find the transition matrix. The series of steps are to:
519  if (anti>0) {
520  // [1] Find the matter Hamiltonian in the mass basis...
521  this->EvalEqn5(twoEH, fU, fUdagg, fdmsqr, E, Gf, Ne);
522 
523  // [2] Find the eigenvalues and sort them.
524  this->EvalEqn22(alpha, beta, gamma, E, Gf, Ne, fdmsqr, fU);
525  this->EvalEqn21(Msqr, alpha, beta, gamma);
526 
527  // Repeat the above, but for vacuum (Ne=0.0) to sort out the order
528  // of the eigenvalues
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);
532 
533  // [3] Evaluate the transition matrix
534  this->EvalEqn11(X, L, E, twoEH, Msqr, dMsqr);
535  this->EvalEqn10(A, fU, X, fUdagg);
536  }
537  else if (anti<0) {
538  // As above, but make required substitutions for anti-neutrinos:
539  // Gf=>-Gf, U=>Ustar, U^dagger=>U^dagger^*=U^T
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);
548  }
549  else abort();
550 
551  // [4] Apply the transition matrix to the matrix we started with...
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];
556  }
557  }
558 }
559 
560 //......................................................................
561 ///
562 /// Do several layers in a row. L and Ne must have the same length
563 ///
564 template <typename T>
565 void _PMNS<T>::PropMatter(const std::list<double>& L,
566  double E,
567  const std::list<double>& Ne,
568  int anti)
569 {
570  if (L.size()!=Ne.size()) abort();
571  auto Li (L.begin());
572  auto Lend(L.end());
573  auto Ni (Ne.begin());
574  for (; Li!=Lend; ++Li, ++Ni) {
575  // For very low densities, use vacumm
576  static const double kRhoCutoff = 1.0E-6; // Ne in moles/cm^3
577  if (*Ni<kRhoCutoff) this->PropVacuum(*Li, E, anti);
578  else this->PropMatter(*Li, E, *Ni, anti);
579  }
580 }
581 
582 //......................................................................
583 ///
584 /// Reset the transition matrix back to the identity matrix from where
585 /// it starts
586 ///
587 template <typename T>
589 {
590  int i, j;
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>();
595  }
596  }
597 }
598 
599 //......................................................................
600 ///
601 /// Compute oscillation probability from flavor i to flavor j
602 ///
603 /// 0 = nue, 1 = numu, 2 = nutau
604 ///
605 template <typename T>
606 T _PMNS<T>::P(int i, int j) const
607 {
608  assert(i>=0 && i<3);
609  assert(j>=0 && j<3);
610  return norm(fA[i][j]);
611 }
612 
613 
614 ////////////////////////////////////////////////////////////////////////
615 
616 // explicit template instantiation for known types
617 
618 template class osc::_PMNS<double>;
619 
620 #ifdef OSCLIB_STAN
621 #include "stan/math/rev/scal.hpp"
622 template class osc::_PMNS<stan::math::var>;
623 #endif
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
void EvalEqn2(complex A[][3], const complex U[][3], const complex Udagg[][3], const T dmsqr[][3], double L, double E)
Definition: PMNS.cxx:241
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
static const double kK2
mole/GeV^2/cm^3 to eV
Definition: PMNSOpt.h:42
double th23
Definition: runWimpSim.h:98
double delta
Definition: runWimpSim.h:98
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
void EvalEqn11(complex X[][3], double L, double E, const complex twoEH[][3], const T Msqr[], const T dMsqr[][3])
Definition: PMNS.cxx:308
void abs(TH1 *hist)
Float_t tmp
Definition: plot.C:36
void PrintDeltaMsqrs() const
Print the matrix of mass-squared differneces.
Definition: PMNS.cxx:172
Double_t beta
double th12
Definition: runWimpSim.h:98
T acos(T number)
Definition: d0nt_math.hpp:54
#define M_PI
Definition: SbMath.h:34
void SortEigenvalues(T dMsqr[][3], const T dmsqr[][3], const T MsqrVac[], T Msqr[])
Definition: PMNS.cxx:441
c2
Definition: demo5.py:33
const double C
static const double kGeV2eV
GeV to eV.
Definition: PMNSOpt.h:43
static constexpr double L
Float_t E
Definition: plot.C:20
void EvalEqn10(complex A[][3], const complex U[][3], const complex X[][3], const complex Udagg[][3])
Definition: PMNS.cxx:292
void EvalEqn22(T &alpha, T &beta, T &gamma, double E, double Gf, double Ne, const T dmsqr[][3], const complex U[][3])
Definition: PMNS.cxx:414
const double a
void Reset()
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
void PropMatter(double L, double E, double Ne, int anti)
Definition: PMNS.cxx:504
Float_t d
Definition: plot.C:236
void PropVacuum(double L, double E, int anti)
Definition: PMNS.cxx:482
const double j
Definition: BetheBloch.cxx:29
void Reset()
Definition: PMNS.cxx:588
void SetMixBWCP(double th1, double th2, double th3, double deltacp)
Definition: PMNS.cxx:139
Oscillation probability calculators.
Definition: Calcs.h:5
OStream cout
Definition: OStream.cxx:6
static const double A
Definition: Units.h:82
T sin(T number)
Definition: d0nt_math.hpp:132
Float_t norm
void One()
const hit & b
Definition: hits.cxx:21
void Multi(complex A[][3], const complex B[][3], const complex C[][3])
Definition: PMNS.cxx:220
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
void SetMix(const T &th12, const T &th23, const T &th13, const T &deltacp)
Definition: PMNS.cxx:97
_PMNS()
Definition: PMNS.cxx:59
c1
Definition: demo5.py:24
void SetDeltaMsqrs(const T &dm21, const T &dm32)
Definition: PMNS.cxx:187
double T
Definition: Xdiff_gwt.C:5
auto one()
Definition: PMNS.cxx:49
void PrintMix() const
Print the oscillation matrix.
Definition: PMNS.cxx:92
T P(int i, int j) const
Definition: PMNS.cxx:606
void EvalEqn5(complex twoEH[][3], const complex U[][3], const complex Udagg[][3], const T dmsqr[][3], double E, double Gf, double Ne)
Definition: PMNS.cxx:265
Float_t X
Definition: plot.C:38
auto zero()
Definition: PMNS.cxx:47
static const double kK1
(1/2)*(1000/hbarc)
Definition: PMNS.cxx:52
double th13
Definition: runWimpSim.h:98
void EvalEqn21(T Msqr[], const T &alpha, const T &beta, const T &gamma)
Definition: PMNS.cxx:368
void DumpMatrix(const complex M[][3]) const
Utility to print a genetrix complex matrix M.
Definition: PMNS.cxx:80