PMNS_Fast.cxx

Go to the documentation of this file.

 
//
// Implementation of oscillations of neutrinos in matter in a
// three-neutrino framework.
//
// jcoelho\@apc.in2p3.fr
 
#include "MatrixDecomp/zheevh3.h"
 
#include "PMNS_Fast.h"
 
using namespace OscProb;
 
//.............................................................................
PMNS_Fast::PMNS_Fast() : PMNS_Base(), fHam() {}
 
//.............................................................................
PMNS_Fast::~PMNS_Fast() {}
 
//.............................................................................
void PMNS_Fast::SetMix(double th12, double th23, double th13, double deltacp)
{
  SetAngle(1, 2, th12);
  SetAngle(1, 3, th13);
  SetAngle(2, 3, th23);
  SetDelta(1, 3, deltacp);
}
 
//.............................................................................
void PMNS_Fast::SetDeltaMsqrs(double dm21, double dm32)
{
  SetDm(2, dm21);
  SetDm(3, dm32 + dm21);
}
 
//.............................................................................
void PMNS_Fast::UpdateHam()
{
  double lv = 2 * kGeV2eV * fEnergy; // 2E in eV
 
  double kr2GNe = kK2 * M_SQRT2 * kGf;
  kr2GNe *= fPath.density * fPath.zoa; // Matter potential in eV
 
  // Finish building Hamiltonian in matter with dimension of eV
  for (int i = 0; i < fNumNus; i++) {
    fHam[i][i] = fHms[i][i] / lv;
    for (int j = i + 1; j < fNumNus; j++) {
      if (!fIsNuBar)
        fHam[i][j] = fHms[i][j] / lv;
      else
        fHam[i][j] = conj(fHms[i][j]) / lv;
    }
  }
  if (!fIsNuBar)
    fHam[0][0] += kr2GNe;
  else
    fHam[0][0] -= kr2GNe;
}
 
//.............................................................................
void PMNS_Fast::SolveHam()
{
  // Do vacuum oscillation in low density
  if (fPath.density < 1.0e-6) {
    SetVacuumEigensystem();
    return;
  }
 
  SolveHamMatter();
}
 
//.............................................................................
void PMNS_Fast::SolveHamMatter()
{
  // Build Hamiltonian
  BuildHms();
 
  // Check if anything changed
  if (fGotES) return;
 
  // Try caching if activated
  if (TryCache()) return;
 
  UpdateHam();
 
  double   fEvalGLoBES[3];
  complexD fEvecGLoBES[3][3];
 
  // Solve Hamiltonian for eigensystem using the GLoBES method
  zheevh3(fHam, fEvecGLoBES, fEvalGLoBES);
 
  // Fill fEval and fEvec vectors from GLoBES arrays
  for (int i = 0; i < fNumNus; i++) {
    fEval[i] = fEvalGLoBES[i];
    for (int j = 0; j < fNumNus; j++) { fEvec[i][j] = fEvecGLoBES[i][j]; }
  }
 
  fGotES = true;
 
  // Fill cache if activated
  FillCache();
}
 
//.............................................................................
void PMNS_Fast::SetVacuumEigensystem()
{
  double   s12, s23, s13, c12, c23, c13;
  complexD idelta(0.0, fDelta[0][2]);
  if (fIsNuBar) idelta = conj(idelta);
 
  s12 = sin(fTheta[0][1]);
  s23 = sin(fTheta[1][2]);
  s13 = sin(fTheta[0][2]);
  c12 = cos(fTheta[0][1]);
  c23 = cos(fTheta[1][2]);
  c13 = cos(fTheta[0][2]);
 
  fEvec[0][0] = c12 * c13;
  fEvec[0][1] = s12 * c13;
  fEvec[0][2] = s13 * exp(-idelta);
 
  fEvec[1][0] = -s12 * c23 - c12 * s23 * s13 * exp(idelta);
  fEvec[1][1] = c12 * c23 - s12 * s23 * s13 * exp(idelta);
  fEvec[1][2] = s23 * c13;
 
  fEvec[2][0] = s12 * s23 - c12 * c23 * s13 * exp(idelta);
  fEvec[2][1] = -c12 * s23 - s12 * c23 * s13 * exp(idelta);
  fEvec[2][2] = c23 * c13;
 
  fEval[0] = 0;
  fEval[1] = fDm[1] / (2 * kGeV2eV * fEnergy);
  fEval[2] = fDm[2] / (2 * kGeV2eV * fEnergy);
}