PMNS_Sterile.cxx

Go to the documentation of this file.

 
 
#include <Eigen/Eigenvalues>
 
#include "PMNS_Sterile.h"
 
using namespace std;
 
using namespace OscProb;
 
//.............................................................................
PMNS_Sterile::PMNS_Sterile(int numNus) : PMNS_Base(numNus), fHam(numNus, numNus)
{
}
 
//.............................................................................
void PMNS_Sterile::UpdateHam()
{
  double rho = fPath.density;
  double zoa = fPath.zoa;
 
  double lv = 2 * kGeV2eV * fEnergy; // 2E in eV
  double kr2GNe =
      kK2 * M_SQRT2 * kGf * rho * zoa; // Electron matter potential in eV
  double kr2GNn = kK2 * M_SQRT2 * kGf * rho * (1 - zoa) /
                  2; // Neutron matter potential in eV
 
  // Finish building Hamiltonian in matter with dimension of eV
  for (int i = 0; i < fNumNus; i++) {
    // Set diagonal elements which are real
    fHam(i, i) = fHms[i][i] / lv;
 
    // Set off-diagonal elements and flip deltaCP if needed
    for (int j = i + 1; j < fNumNus; j++) {
      // Obs: fHam is lower tirangular while fHms is upper triangular
      if (!fIsNuBar)
        fHam(j, i) = conj(fHms[i][j]) / lv;
      else
        fHam(j, i) = fHms[i][j] / lv;
    }
 
    // Subtract NC coherent forward scattering from sterile neutrinos.
    // See arXiv:hep-ph/0606054v3, eq. 3.15, for example.
    if (i > 2) {
      if (!fIsNuBar)
        fHam(i, i) += kr2GNn;
      else
        fHam(i, i) -= kr2GNn;
    }
  }
 
  // Add nue CC coherent forward scattering.
  if (!fIsNuBar)
    fHam(0, 0) += kr2GNe;
  else
    fHam(0, 0) -= kr2GNe;
}
 
//.............................................................................
template <typename T> void PMNS_Sterile::SolveEigenSystem()
{
  Eigen::Ref<T> H(fHam);
 
  Eigen::SelfAdjointEigenSolver<T> eigensolver(H);
 
  // Fill fEval and fEvec vectors from GSL objects
  for (int i = 0; i < fNumNus; i++) {
    fEval[i] = eigensolver.eigenvalues()(i);
    for (int j = 0; j < fNumNus; j++) {
      fEvec[i][j] = eigensolver.eigenvectors()(i, j);
    }
  }
}
 
//.............................................................................
void PMNS_Sterile::SolveHam()
{
  // Build Hamiltonian
  BuildHms();
 
  // Check if anything changed
  if (fGotES) return;
 
  // Try caching if activated
  if (TryCache()) return;
 
  UpdateHam();
 
  // Choose the appropriate MatrixType for the solver
  if (fNumNus == 4)
    SolveEigenSystem<Eigen::Matrix4cd>();
  else if (fNumNus == 3)
    SolveEigenSystem<Eigen::Matrix3cd>();
  else if (fNumNus == 2)
    SolveEigenSystem<Eigen::Matrix2cd>();
  else
    SolveEigenSystem<Eigen::MatrixXcd>();
 
  // Mark eigensystem as solved
  fGotES = true;
 
  // Fill cache if activated
  FillCache();
}