#include <PMNS.h>

Public Member Functions
	PMNS ()

	PMNS (double th12, double th23, double th13, double deltacp, double dms12, double dms23)

void	PrintMix () const

void	PrintDeltaMsqrs () const

double	P (int i, int j) const

void	SetMix (double th12, double th23, double th13, double deltacp)

void	SetMixBWCP (double th1, double th2, double th3, double deltacp)

void	SetDeltaMsqrs (double dm21, double dm32)

void	PropVacuum (double L, double E, int anti)

void	PropMatter (double L, double E, double Ne, int anti)

void	PropMatter (const std::list< double > &L, double E, const std::list< double > &Ne, int anti)

void	Reset ()

Private Types
typedef std::complex< double >	complex

Private Member Functions
void	Multi (complex A[][3], const complex B[][3], const complex C[][3])

void	EvalEqn2 (complex A[][3], const complex U[][3], const complex Udagg[][3], const double dmsqr[][3], double L, double E)

void	EvalEqn5 (complex twoEH[][3], const complex U[][3], const complex Udagg[][3], const double dmsqr[][3], double E, double Gf, double Ne)

void	EvalEqn10 (complex A[][3], const complex U[][3], const complex X[][3], const complex Udagg[][3])

void	EvalEqn11 (complex X[][3], double L, double E, const complex twoEH[][3], const double Msqr[], const double dMsqr[][3])

void	EvalEqn21 (double Msqr[], double alpha, double beta, double gamma)

void	EvalEqn22 (double &alpha, double &beta, double &gamma, double E, double Gf, double Ne, const double dmsqr[][3], const complex U[][3])

void	SortEigenvalues (double dMsqr[][3], const double dmsqr[][3], const double MsqrVac[], double Msqr[])

void	DumpMatrix (const complex M[][3]) const

Private Attributes
double	fdmsqr [3][3]

complex	fU [3][3]

complex	fUdagg [3][3]

complex	fUstar [3][3]

complex	fUtran [3][3]

complex	fA [3][3]

Detailed Description

Definition at line 38 of file PMNS.h.

Member Typedef Documentation

typedef std::complex<double> osc::PMNS::complex

private

Definition at line 73 of file PMNS.h.

Constructor & Destructor Documentation

PMNS::PMNS ( )

Definition at line 52 of file PMNS.cxx.

 {
   this->SetMix(0.,0.,0.,0.);
   this->SetDeltaMsqrs(0.,0.);
   this->Reset();
 }

PMNS::PMNS	(	double	th12,
		double	th23,
		double	th13,
		double	deltacp,
		double	dms12,
		double	dms23
	)

Definition at line 61 of file PMNS.cxx.

 {
   this->SetMix(th12, th23, th13, deltacp);
   this->SetDeltaMsqrs(dms21, dms32);
   this->Reset();
 }

Member Function Documentation

void PMNS::DumpMatrix ( const complex M[][3] ) const

private

Definition at line 71 of file PMNS.cxx.

 {
   std::cout 
     <<"| "<<M[0][0]<<"\t"<<M[0][1]<<"\t"<<M[0][2]<<" |\n"
     <<"| "<<M[1][0]<<"\t"<<M[1][1]<<"\t"<<M[1][2]<<" |\n"
     <<"| "<<M[2][0]<<"\t"<<M[2][1]<<"\t"<<M[2][2]<<" |\n"
     <<std::endl;
 }

void PMNS::EvalEqn10	(	complex	A[][3],
		const complex	U[][3],
		const complex	X[][3],
		const complex	Udagg[][3]
	)

private

Definition at line 242 of file PMNS.cxx.

 {
   complex tmp[3][3];
   this->Multi(tmp, X, Udagg);
   this->Multi(A,   U, tmp);
 }

void PMNS::EvalEqn11	(	complex	X[][3],
		double	L,
		double	E,
		const complex	twoEH[][3],
		const double	Msqr[],
		const double	dMsqr[][3]
	)

private

Definition at line 253 of file PMNS.cxx.

 {
   // The identity matrix
   static const double One[3][3] = {{1.,0.,0.},
                                    {0.,1.,0.},
                                    {0.,0.,1.}
   };
   register int a, b, k;
   complex phase;
   complex EHM0[3][3];
   complex EHM1[3][3];
   complex EHM2[3][3];
   complex EHM21[3][3];
   complex EHM20[3][3];
   complex EHM10[3][3];
 
   // There are three matrices which can apper inside the product on
   // j!=k. Calculate them before hand
   for (a=0; a<3; ++a) {
     for (b=0; b<3; ++b) {
       EHM0[a][b] = twoEH[a][b]-Msqr[0]*One[a][b];
       EHM1[a][b] = twoEH[a][b]-Msqr[1]*One[a][b];
       EHM2[a][b] = twoEH[a][b]-Msqr[2]*One[a][b];
     }
   }
   this->Multi(EHM21,EHM2,EHM1);
   this->Multi(EHM20,EHM2,EHM0);
   this->Multi(EHM10,EHM1,EHM0);
 
   // Refer to Msqr_j as dMsqr[j][0] since only mass differences matter
   for (a=0; a<3; ++a) {
     for (b=0; b<3; ++b) {
       X[a][b] = zero;
       for (k=0; k<3; ++k) {
         phase = exp(complex(0.0,-kK1*dMsqr[k][0]*L/E));
         if (k==0) {
           X[a][b] += (EHM21[a][b]/(dMsqr[k][2]*dMsqr[k][1]))*phase;
         }
         else if (k==1) { 
           X[a][b] += (EHM20[a][b]/(dMsqr[k][2]*dMsqr[k][0]))*phase;
         }
         else if (k==2) {
           X[a][b] += (EHM10[a][b]/(dMsqr[k][1]*dMsqr[k][0]))*phase;
         } // Switch for product on j!=k
       } // Sum on k
     } // Loop on b
   } // Loop on a
 }

void PMNS::EvalEqn2	(	complex	A[][3],
		const complex	U[][3],
		const complex	Udagg[][3],
		const double	dmsqr[][3],
		double	L,
		double	E
	)

private

Definition at line 202 of file PMNS.cxx.

 {
   register int a, b, i;
   for (a=0; a<3; ++a) {
     for (b=0; b<3; ++b) {
       A[a][b] = zero;
       for (i=0; i<3; ++i) {
         complex phase(0.0,-kK1*dmsqr[i][0]*L/E);
         A[a][b] += U[a][i]*exp(phase)*Udagg[i][b];
       }
     }
   }
 }

void PMNS::EvalEqn21	(	double	Msqr[],
		double	alpha,
		double	beta,
		double	gamma
	)

private

Definition at line 312 of file PMNS.cxx.

 {
   double arg;      // Argument of the acos()
   double theta0;   // First of the three roots of acos()
   double theta1;   // Second of the three roots of acos()
   double theta2;   // Third of the three roots of acos()
   double fac;      // Factor in front of cos() terms
   double constant; // Offset for all eigenvalues
   static const double k2PiOver3 = 2.0*M_PI/3.0;
   double alpha2           = alpha*alpha;
   double alpha3           = alpha*alpha2;
   double alpha2Minus3beta = alpha2-3.0*beta;
 
   arg =
     (2.0*alpha3 - 9.0*alpha*beta + 27.0*gamma)/
     (2.0*pow(alpha2Minus3beta,1.5));
 
   // Occasionally round off errors mean that arg wanders outside of
   // its allowed range. If its way off (1 part per billion), stop the
   // program. Otherwise, set it to its real value.
   if (fabs(arg)>1.0) {
     if (fabs(arg-1.0)>1.E-9) abort();
     if (arg<-1.0) arg = -1.00;
     if (arg>+1.0) arg = +1.00;
   }
   
   // The three roots, find the first by computing the acos() the
   // others are nearby
   theta0 = acos(arg)/3.0;
   theta1 = theta0-k2PiOver3;
   theta2 = theta0+k2PiOver3;
   
   // The multiplier and offset
   fac      = -2.0*sqrt(alpha2Minus3beta)/3.0;
   constant = -alpha/3.0; // The constant offset m1^2 is irrelevant
   
   // And the eigenvalues themselves
   Msqr[0] = fac*cos(theta0) + constant;
   Msqr[1] = fac*cos(theta1) + constant;
   Msqr[2] = fac*cos(theta2) + constant;
 }

void PMNS::EvalEqn22	(	double &	alpha,
		double &	beta,
		double &	gamma,
		double	E,
		double	Gf,
		double	Ne,
		const double	dmsqr[][3],
		const complex	U[][3]
	)

private

Definition at line 357 of file PMNS.cxx.

 {
   // 2*sqrt(2)*Gf*Ne*E in units of eV^2
   double k2r2EGNe = kK2*2.0*M_SQRT2*Gf*Ne*(kGeV2eV*E);
   
   alpha = k2r2EGNe + dmsqr[0][1] + dmsqr[0][2];
   
   beta  =
     dmsqr[0][1]*dmsqr[0][2] + 
     k2r2EGNe*(dmsqr[0][1]*(1.0-norm(U[0][1])) + 
               dmsqr[0][2]*(1.0-norm(U[0][2])));
 
   gamma = k2r2EGNe*dmsqr[0][1]*dmsqr[0][2]*norm(U[0][0]);
 }

void PMNS::EvalEqn5	(	complex	twoEH[][3],
		const complex	U[][3],
		const complex	Udagg[][3],
		const double	dmsqr[][3],
		double	E,
		double	Gf,
		double	Ne
	)

private

Definition at line 222 of file PMNS.cxx.

 {
   register int j, k;
   double       k2r2GNeE = kK2*2.0*M_SQRT2*Gf*Ne*(kGeV2eV*E);
   for (k=0; k<3; ++k) {
     for (j=0; j<3; ++j) {
       twoEH[k][j] = zero;
       if (k==j) twoEH[k][j] = dmsqr[j][0];
       twoEH[k][j] -= k2r2GNeE*U[0][j]*Udagg[k][0];
     }
   }
 }

void PMNS::Multi	(	complex	A[][3],
		const complex	B[][3],
		const complex	C[][3]
	)

private

Definition at line 188 of file PMNS.cxx.

 {
   register int i, j, k;
   for (i=0; i<3; ++i) {
     for (j=0; j<3; ++j) {
       A[i][j] = zero;
       for (k=0; k<3; ++k) {
         A[i][j] += B[i][k]*C[k][j];
       }
     }
   }
 }

double PMNS::P	(	int	i,
		int	j
	)		const

Definition at line 519 of file PMNS.cxx.

 {
   assert(i>=0 && i<3);
   assert(j>=0 && j<3);
   return norm(fA[i][j]);
 }

void PMNS::PrintDeltaMsqrs ( ) const

Definition at line 149 of file PMNS.cxx.

 {
   std::cout 
     <<"|"<<fdmsqr[0][0]<<"\t"<<fdmsqr[0][1]<<"\t"<<fdmsqr[0][2]<<"|\n" 
     <<"|"<<fdmsqr[1][0]<<"\t"<<fdmsqr[1][1]<<"\t"<<fdmsqr[1][2]<<"|\n" 
     <<"|"<<fdmsqr[2][0]<<"\t"<<fdmsqr[2][1]<<"\t"<<fdmsqr[2][2]<<"|"
     << std::endl;
 }

void PMNS::PrintMix ( ) const

Definition at line 82 of file PMNS.cxx.

82 { this->DumpMatrix(fU); }

osc::PMNS::fU

complex fU[3][3]

Definition: PMNS.h:127

osc::PMNS::DumpMatrix

void DumpMatrix(const complex M[][3]) const

Definition: PMNS.cxx:71

void PMNS::PropMatter	(	double	L,
		double	E,
		double	Ne,
		int	anti
	)

Definition at line 432 of file PMNS.cxx.

 {
   static const double  Gf = 1.166371E-5; // G_F in units of GeV^-2
   register int i, j;
   complex twoEH[3][3];
   complex X[3][3];
   double  Msqr[3];
   double  MsqrV[3];
   double  dMsqr[3][3];
   double  alpha,  beta,  gamma;
   double  alpha0, beta0, gamma0;
   complex A[3][3];
   complex Aout[3][3];
   
   // Find the transition matrix. The series of steps are to:
   if (anti>0) {
     // [1] Find the matter Hamiltonian in the mass basis...
     this->EvalEqn5(twoEH, fU, fUdagg, fdmsqr, E, Gf, Ne);
 
     // [2] Find the eigenvalues and sort them.
     this->EvalEqn22(alpha, beta, gamma, E, Gf, Ne, fdmsqr, fU);
     this->EvalEqn21(Msqr,  alpha, beta, gamma);
 
     // Repeat the above, but for vacuum (Ne=0.0) to sort out the order
     // of the eigenvalues
     this->EvalEqn22(alpha0, beta0, gamma0, E, 0.0, 0.0, fdmsqr, fU);
     this->EvalEqn21(MsqrV, alpha0, beta0, gamma0);
     this->SortEigenvalues(dMsqr, fdmsqr, MsqrV, Msqr);
 
     // [3] Evaluate the transition matrix
     this->EvalEqn11(X, L, E, twoEH, Msqr, dMsqr);
     this->EvalEqn10(A, fU, X, fUdagg);
   }
   else if (anti<0) {
     // As above, but make required substitutions for anti-neutrinos:
     // Gf=>-Gf, U=>Ustar, U^dagger=>U^dagger^*=U^T
     this->EvalEqn5(twoEH, fUstar, fUtran, fdmsqr, E, -Gf, Ne);
     this->EvalEqn22(alpha, beta, gamma, E, -Gf, Ne, fdmsqr, fUstar);
     this->EvalEqn21(Msqr,  alpha, beta, gamma);
     this->EvalEqn22(alpha0, beta0, gamma0, E, 0.0, 0.0, fdmsqr, fUstar);
     this->EvalEqn21(MsqrV, alpha0, beta0, gamma0);
     this->SortEigenvalues(dMsqr, fdmsqr, MsqrV, Msqr);
     this->EvalEqn11(X, L, E, twoEH, Msqr, dMsqr);
     this->EvalEqn10(A, fUstar, X, fUtran);
   }
   else abort();
   
   // [4] Apply the transition matrix to the matrix we started with...
   this->Multi(Aout, A, fA);
   for (i=0; i<3; ++i) {
     for (j=0; j<3; ++j) {
       fA[i][j] = Aout[i][j];
     }
   }
 }

void PMNS::PropMatter	(	const std::list< double > &	L,
		double	E,
		const std::list< double > &	Ne,
		int	anti
	)

Definition at line 489 of file PMNS.cxx.

 {
   if (L.size()!=Ne.size()) abort();
   std::list<double>::const_iterator Li  (L.begin());
   std::list<double>::const_iterator Lend(L.end());
   std::list<double>::const_iterator Ni  (Ne.begin());
   for (; Li!=Lend; ++Li, ++Ni) {
     // For very low densities, use vacumm
     static const double kRhoCutoff = 1.0E-6; // Ne in moles/cm^3
     if (*Ni<kRhoCutoff) this->PropVacuum(*Li, E, anti);
     else                this->PropMatter(*Li, E, *Ni, anti);
   }
 }

void PMNS::PropVacuum	(	double	L,
		double	E,
		int	anti
	)

Definition at line 415 of file PMNS.cxx.

 {
   register int i, j;
   complex A[3][3];
   complex Aout[3][3];
 
   if      (anti>0) this->EvalEqn2(A, fU,     fUdagg, fdmsqr, L, E);
   else if (anti<0) this->EvalEqn2(A, fUstar, fUtran, fdmsqr, L, E);
   else abort();
   this->Multi(Aout, A, fA);
   for (i=0; i<3; ++i) {
     for (j=0; j<3; ++j) {
       fA[i][j] = Aout[i][j];
     }
   }
 }

void PMNS::Reset ( )

Definition at line 507 of file PMNS.cxx.

 {
   register int i, j;
   for (i=0; i<3; ++i) {
     for (j=0; j<3; ++j) {
       if (i==j) fA[i][j] = one;
       else      fA[i][j] = zero;
     }
   }
 }

void PMNS::SetDeltaMsqrs	(	double	dm21,
		double	dm32
	)

Definition at line 159 of file PMNS.cxx.

 {
   double eps = 5.0E-9;
   double msqr[3];
   
   msqr[0] = 0.0;
   msqr[1] = dm21;
   msqr[2] = dm21+dm32;
   
   // Degeneracies cause problems with diagonalization, so break them
   // ever so slightly
   if (dm21==0.0) {msqr[0] -= 0.5*eps; msqr[1] += 0.5*eps; }
   if (dm32==0.0) {msqr[1] -= 0.5*eps; msqr[2] += 0.5*eps; }
 
   // Assign the mass splittings fdmsqr_ij = msqr_i - msqr_j by
   // convention
   for (register int i=0; i<3; ++i) {
     for (register int j=0; j<3; ++j) {
       // A somewhat subtle point: Barger et al. refer to the sign of
       // m1-m2 being positive which implies dm^2_12>0 and
       // dm^2_21<0. The labeling in more common use is to assume m3 is
       // heaviest such that dm_12<0 and dm_21>0. Rather than reverse
       // all the indices in all the equations, flip the signs here.
       fdmsqr[i][j] = -(msqr[i] - msqr[j]);
     }
   }
 }

void PMNS::SetMix	(	double	th12,
		double	th23,
		double	th13,
		double	deltacp
	)

Definition at line 86 of file PMNS.cxx.

 {
   register int i, j;
   double  s12, s23, s13, c12, c23, c13;
   complex idelta(0.0,deltacp);
   
   s12 = sin(th12);  s23 = sin(th23);  s13 = sin(th13);
   c12 = cos(th12);  c23 = cos(th23);  c13 = cos(th13);
   
   fU[0][0] =  c12*c13;
   fU[0][1] =  s12*c13;
   fU[0][2] =  s13*exp(-idelta);
   
   fU[1][0] = -s12*c23-c12*s23*s13*exp(idelta);
   fU[1][1] =  c12*c23-s12*s23*s13*exp(idelta);
   fU[1][2] =  s23*c13;
   
   fU[2][0] =  s12*s23-c12*c23*s13*exp(idelta);
   fU[2][1] = -c12*s23-s12*c23*s13*exp(idelta);
   fU[2][2] =  c23*c13;
   
   // Compute derived forms of the mixing matrix
   for (i=0; i<3; ++i) {
     for (j=0; j<3; ++j) {
       fUtran[i][j] = fU[j][i];
       fUstar[i][j] = conj(fU[i][j]);
       fUdagg[i][j] = conj(fU[j][i]);
     }
   }
 }

void PMNS::SetMixBWCP	(	double	th1,
		double	th2,
		double	th3,
		double	deltacp
	)

Definition at line 117 of file PMNS.cxx.

 {
   int i, j;
   double s1, s2, s3, c1, c2, c3;
   complex id(0.0,d);
   s1 = sin(th1);  s2 = sin(th2);  s3 = sin(th3);
   c1 = cos(th1);  c2 = cos(th2);  c3 = cos(th3);
   
   fU[0][0] =  c1;
   fU[0][1] =  s1*c3;
   fU[0][2] =  s1*s3;
 
   fU[1][0] = -s1*c2;
   fU[1][1] =  c1*c2*c3+s2*s3*exp(id);
   fU[1][2] =  c1*c2*s3-s2*c3*exp(id);
 
   fU[2][0] = -s1*s2;
   fU[2][1] =  c1*s2*c3-c2*s3*exp(id);
   fU[2][2] =  c1*s2*s3+c2*c3*exp(id);
 
   // Compute derived forms of the mixing matrix
   for (i=0; i<3; ++i) {
     for (j=0; j<3; ++j) {
       fUtran[i][j] = fU[j][i];
       fUstar[i][j] = conj(fU[i][j]);
       fUdagg[i][j] = conj(fU[j][i]);
     }
   }
 }

void PMNS::SortEigenvalues	(	double	dMsqr[][3],
		const double	dmsqr[][3],
		const double	MsqrVac[],
		double	Msqr[]
	)

private

Definition at line 380 of file PMNS.cxx.

 {
   int i, j, k;
   double best, delta;
   double MsqrTmp[3] = {-99.9E9,-99.9E9,-99.9E9};
   int flg[3] = {0,0,0};
 
   // Attempt to figure out which of the eigenvalues match between
   // dmsqr and MsqrVac
   for (i=0; i<3; ++i) {
     best =  1.E30;
     k    = -1;
     for (j=0; j<3; ++j) {
       delta = fabs(MsqrVac[i] - dmsqr[j][0]);
       if (delta<best) { best = delta; k = j; }
     }
     if (best>1.E-9) abort();
     if (k<0 || k>2) abort();
     flg[k] = 1;
     MsqrTmp[i] = Msqr[k];
   }
   // Check that each eigenvalue is used
   for (i=0; i<3; ++i) if (flg[i]!=1) abort();
   
   for (i=0; i<3; ++i) {
     Msqr[i] = MsqrTmp[i];
     for (j=0; j<3; ++j) {
       dMsqr[i][j] = (MsqrTmp[i] - MsqrTmp[j]);
     }
   }
 }

Member Data Documentation

complex osc::PMNS::fA[3][3]

private

Definition at line 131 of file PMNS.h.

double osc::PMNS::fdmsqr[3][3]

private

Definition at line 126 of file PMNS.h.

complex osc::PMNS::fU[3][3]

private

Definition at line 127 of file PMNS.h.

complex osc::PMNS::fUdagg[3][3]

private

Definition at line 128 of file PMNS.h.

complex osc::PMNS::fUstar[3][3]

private

Definition at line 129 of file PMNS.h.

complex osc::PMNS::fUtran[3][3]

private

Definition at line 130 of file PMNS.h.

The documentation for this class was generated from the following files:

g4lbne/ProductionScripts/OscLib/PMNS.h
g4lbne/ProductionScripts/OscLib/PMNS.cxx

Public Member Functions

Private Types

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation