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