PMNS_OQS.cxx

Go to the documentation of this file.

 
 
#include <unsupported/Eigen/MatrixFunctions>
 
#include "PMNS_OQS.h"
#include "exceptions.h"
 
#include <iostream>
 
using namespace std;
using namespace OscProb;
 
constexpr int SU3_DIM = PMNS_OQS::SU3_DIM;
 
//.............................................................................
PMNS_OQS::PMNS_OQS()
    : PMNS_DensityMatrix(), fUM(3, vectorC(3, 0)), fHVMB(3, vectorC(3, 0)),
      fHGM(SU3_DIM, vectorD(SU3_DIM, 0)), fD(SU3_DIM, vectorD(SU3_DIM, 0)),
      fcos(SU3_DIM, vectorD(SU3_DIM, 1)), fa(SU3_DIM, 0), fR(SU3_DIM),
      fRt(SU3_DIM), fM(8, 8)
{
  SetPower(0);
  fBuiltDissipator = false;
}
 
//.............................................................................
PMNS_OQS::~PMNS_OQS() {}
 
//.............................................................................
void PMNS_OQS::SetPower(int n) { fPower = n; }
 
//.............................................................................
void PMNS_OQS::SetDecoElement(int i, double val)
{
  THROW_ON_INVALID_ARG(i > 0, i);
  THROW_ON_INVALID_ARG(i < SU3_DIM, i, SU3_DIM);
 
  fBuiltDissipator *= (fa[i] == abs(val));
  fa[i] = abs(val);
}
 
//.............................................................................
void PMNS_OQS::SetDecoAngle(int i, int j, double th)
{
  THROW_ON_INVALID_ARG(i > 0, i);
  THROW_ON_INVALID_ARG(i < SU3_DIM, i, SU3_DIM);
  THROW_ON_INVALID_ARG(j > 0, i);
  THROW_ON_INVALID_ARG(j < SU3_DIM, j, SU3_DIM);
  THROW_ON_INVALID_ARG(i != j, i, j);
 
  double val = cos(th);
  fBuiltDissipator *= (fcos[i][j] == val);
  fcos[i][j] = val;
  fcos[j][i] = val;
}
 
//.............................................................................
int PMNS_OQS::GetPower() { return fPower; }
 
//.............................................................................
double PMNS_OQS::GetDecoElement(int i)
{
  THROW_ON_INVALID_ARG(i > 0, i);
  THROW_ON_INVALID_ARG(i < SU3_DIM, i, SU3_DIM);
 
  return fa[i];
}
 
//.............................................................................
double PMNS_OQS::GetDecoAngle(int i, int j)
{
  THROW_ON_INVALID_ARG(i > 0, i);
  THROW_ON_INVALID_ARG(i < SU3_DIM, i, SU3_DIM);
  THROW_ON_INVALID_ARG(j > 0, i);
  THROW_ON_INVALID_ARG(j < SU3_DIM, j, SU3_DIM);
  THROW_ON_INVALID_ARG(i != j, i, j);
 
  return acos(fcos[i][j]);
}
 
//.............................................................................
double PMNS_OQS::GetHGM(int i, int j)
{
  THROW_ON_INVALID_ARG(i >= 0, i);
  THROW_ON_INVALID_ARG(i < SU3_DIM, i, SU3_DIM);
  THROW_ON_INVALID_ARG(j >= 0, i);
  THROW_ON_INVALID_ARG(j < SU3_DIM, j, SU3_DIM);
 
  return fHGM[i][j];
}
 
//.............................................................................
double PMNS_OQS::GetDissipatorElement(int i, int j)
{
  THROW_ON_INVALID_ARG(i >= 0, i);
  THROW_ON_INVALID_ARG(i < SU3_DIM, i, SU3_DIM);
  THROW_ON_INVALID_ARG(j >= 0, i);
  THROW_ON_INVALID_ARG(j < SU3_DIM, j, SU3_DIM);
 
  return fD[i][j];
}
 
//.............................................................................
void PMNS_OQS::SetIsNuBar(bool isNuBar)
{
  fBuiltHms *= (fIsNuBar == isNuBar);
  PMNS_Base::SetIsNuBar(isNuBar);
}
 
//.............................................................................
void PMNS_OQS::BuildHms()
{
  if (fBuiltHms) return;
  SetVacuumEigensystem();
  for (int i = 0; i < 3; i++) {
    for (int j = 0; j < 3; j++) { fUM[i][j] = conj(fEvec[i][j]); }
  }
  PMNS_Base::BuildHms();
}
 
//.............................................................................
void PMNS_OQS::BuildHVMB(NuPath p)
{
  // Set the neutrino path
  SetCurPath(p);
 
  double kr2GNe = kK2 * M_SQRT2 * kGf;
  kr2GNe *= fPath.density * fPath.zoa; // Matter potential in eV
 
  double Ve = 0;
 
  if (!fIsNuBar)
    Ve = kr2GNe;
  else
    Ve = -kr2GNe;
 
  BuildHms();
 
  for (int i = 0; i < 3; i++) {
    for (int j = i; j < 3; j++) {
      fHVMB[i][j] = conj(fUM[0][i]) * Ve * fUM[0][j];
      if (i > j) fHVMB[j][i] = conj(fHVMB[i][j]);
    }
  }
 
  double lv = 2 * kGeV2eV * fEnergy; // 2E in eV
 
  // add mass terms
  fHVMB[1][1] += fDm[1] / lv;
  fHVMB[2][2] += fDm[2] / lv;
}
 
//.............................................................................
void get_GM(const matrixC& A, vectorD& v)
{
  v[0] = real(A[0][0] + A[1][1] + A[2][2]) / sqrt(6);
  v[3] = real(A[0][0] - A[1][1]) / 2;
  v[8] = real(A[0][0] + A[1][1] - 2. * A[2][2]) / sqrt(12);
 
  v[1] = real(A[0][1]);
  v[4] = real(A[0][2]);
  v[6] = real(A[1][2]);
 
  v[2] = -imag(A[0][1]);
  v[5] = -imag(A[0][2]);
  v[7] = -imag(A[1][2]);
}
 
//.............................................................................
void get_SU3(const vectorD& v, matrixC& A)
{
  A[0][0] = v[0] * sqrt(2 / 3.) + v[3] + v[8] / sqrt(3);
  A[1][1] = v[0] * sqrt(2 / 3.) - v[3] + v[8] / sqrt(3);
  A[2][2] = v[0] * sqrt(2 / 3.) - 2 * v[8] / sqrt(3);
 
  A[1][0] = complexD(v[1], v[2]);
  A[2][0] = complexD(v[4], v[5]);
  A[2][1] = complexD(v[6], v[7]);
 
  for (int i = 0; i < 3; i++) {
    for (int j = 0; j < i; j++) { A[j][i] = conj(A[i][j]); }
  }
}
 
//.............................................................................
void get_GMOP(const matrixC& A, matrixD& B)
{
  B[1][2] = real(A[0][0] - A[1][1]);
  B[1][3] = 2 * imag(A[0][1]);
  B[1][4] = -imag(A[1][2]);
  B[1][5] = -real(A[1][2]);
  B[1][6] = -imag(A[0][2]);
  B[1][7] = -real(A[0][2]);
 
  B[2][3] = 2 * real(A[0][1]);
  B[2][4] = real(A[1][2]);
  B[2][5] = -imag(A[1][2]);
  B[2][6] = -real(A[0][2]);
  B[2][7] = imag(A[0][2]);
 
  B[3][4] = -imag(A[0][2]);
  B[3][5] = -real(A[0][2]);
  B[3][6] = imag(A[1][2]);
  B[3][7] = real(A[1][2]);
 
  B[4][5] = real(A[0][0] - A[2][2]);
  B[4][6] = -imag(A[0][1]);
  B[4][7] = real(A[0][1]);
  B[4][8] = sqrt(3) * imag(A[0][2]);
 
  B[5][6] = -real(A[0][1]);
  B[5][7] = -imag(A[0][1]);
  B[5][8] = sqrt(3) * real(A[0][2]);
 
  B[6][7] = real(A[1][1] - A[2][2]);
  B[6][8] = sqrt(3) * imag(A[1][2]);
 
  B[7][8] = sqrt(3) * real(A[1][2]);
 
  for (int i = 1; i < SU3_DIM; ++i) {
    for (int j = i + 1; j < SU3_DIM; ++j) { B[j][i] = -B[i][j]; }
  }
}
 
//.............................................................................
void PMNS_OQS::BuildHGM(NuPath p)
{
  BuildHVMB(p);
  get_GMOP(fHVMB, fHGM);
}
 
//.............................................................................
void PMNS_OQS::BuildDissipator()
{
  if (fBuiltDissipator) return;
 
  double aa[SU3_DIM][SU3_DIM];
  for (int i = 1; i < SU3_DIM; i++) {
    for (int j = i; j < SU3_DIM; j++) {
      aa[i][j] = fa[i] * fa[j];
      if (i == 8) aa[i][j] *= sqrt(3);
      if (j == 8) aa[i][j] *= sqrt(3);
      if (i < j) aa[i][j] *= fcos[i][j];
    }
  }
  double sum12   = aa[1][1] + aa[2][2];
  double sum45   = aa[4][4] + aa[5][5];
  double sum67   = aa[6][6] + aa[7][7];
  double gamma21 = aa[3][3];
  double gamma31 = (aa[3][3] + aa[8][8] + 2 * aa[3][8]) / 4;
  double gamma32 = (aa[3][3] + aa[8][8] - 2 * aa[3][8]) / 4;
 
  fD[1][1] = gamma21 + aa[2][2] + (sum45 + sum67) / 4;
  fD[2][2] = gamma21 + aa[1][1] + (sum45 + sum67) / 4;
  fD[3][3] = sum12 + (sum45 + sum67) / 4;
 
  fD[4][4] = gamma31 + aa[5][5] + (sum12 + sum67) / 4;
  fD[5][5] = gamma31 + aa[4][4] + (sum12 + sum67) / 4;
  fD[6][6] = gamma32 + aa[7][7] + (sum12 + sum45) / 4;
  fD[7][7] = gamma32 + aa[6][6] + (sum12 + sum45) / 4;
 
  fD[8][8] = (sum45 + sum67) * 3 / 4;
  fD[3][8] = (sum45 - sum67) * sqrt(3) / 4;
 
  fD[1][2] = -aa[1][2];
  fD[1][3] = -aa[1][3];
  fD[2][3] = -aa[2][3];
  fD[4][5] = -aa[4][5];
  fD[6][7] = -aa[6][7];
 
  fD[1][8] = (aa[4][6] + aa[5][7]) * sqrt(3) / 2;
  fD[2][8] = (aa[5][6] - aa[4][7]) * sqrt(3) / 2;
 
  fD[4][6] = -(aa[4][6] - 2 * aa[1][8] - 3 * aa[5][7]) / 4;
  fD[4][7] = -(aa[4][7] + 2 * aa[2][8] + 3 * aa[5][6]) / 4;
  fD[5][6] = -(aa[5][6] - 2 * aa[2][8] + 3 * aa[4][7]) / 4;
  fD[5][7] = -(aa[5][7] - 2 * aa[1][8] - 3 * aa[4][6]) / 4;
 
  fD[1][4] = -(aa[1][4] - 3 * (aa[2][5] - aa[3][6]) + aa[6][8]) / 4;
  fD[1][5] = -(aa[1][5] + 3 * (aa[2][4] + aa[3][7]) + aa[7][8]) / 4;
  fD[1][6] = -(aa[1][6] + 3 * (aa[2][7] - aa[3][4]) + aa[4][8]) / 4;
  fD[1][7] = -(aa[1][7] - 3 * (aa[2][6] + aa[3][5]) + aa[5][8]) / 4;
  fD[2][4] = -(aa[2][4] + 3 * (aa[1][5] - aa[3][7]) - aa[7][8]) / 4;
  fD[2][5] = -(aa[2][5] - 3 * (aa[1][4] - aa[3][6]) + aa[6][8]) / 4;
  fD[2][6] = -(aa[2][6] - 3 * (aa[1][7] + aa[3][5]) + aa[5][8]) / 4;
  fD[2][7] = -(aa[2][7] + 3 * (aa[1][6] + aa[3][4]) - aa[4][8]) / 4;
  fD[3][4] = -(aa[3][4] - 3 * (aa[1][6] - aa[2][7]) + aa[4][8]) / 4;
  fD[3][5] = -(aa[3][5] - 3 * (aa[1][7] + aa[2][6]) + aa[5][8]) / 4;
  fD[3][6] = -(aa[3][6] + 3 * (aa[1][4] + aa[2][5]) - aa[6][8]) / 4;
  fD[3][7] = -(aa[3][7] + 3 * (aa[1][5] - aa[2][4]) - aa[7][8]) / 4;
 
  fD[4][8] = -(aa[1][6] - aa[2][7] + aa[3][4] + aa[4][8]) * sqrt(3) / 4;
  fD[5][8] = -(aa[1][7] + aa[2][6] + aa[3][5] + aa[5][8]) * sqrt(3) / 4;
  fD[6][8] = -(aa[1][4] + aa[2][5] - aa[3][6] + aa[6][8]) * sqrt(3) / 4;
  fD[7][8] = -(aa[1][5] - aa[2][4] - aa[3][7] + aa[7][8]) * sqrt(3) / 4;
 
  for (int j = 1; j < SU3_DIM; j++) {
    for (int k = j; k < SU3_DIM; k++) {
      fD[j][k] = -fD[j][k] * kGeV2eV;
      fD[k][j] = fD[j][k];
    }
  }
 
  fBuiltDissipator = true;
}
 
//.............................................................................
void PMNS_OQS::BuildM(NuPath p)
{
  BuildHGM(p);
  BuildDissipator();
  double energyCorr = pow(fEnergy, fPower);
  for (int k = 1; k < SU3_DIM; ++k) {
    for (int j = 1; j < SU3_DIM; ++j) {
      fM(k - 1, j - 1) = fHGM[k][j] + fD[k][j] * energyCorr;
    }
  }
}
 
//.............................................................................
void PMNS_OQS::BuildR() { get_GM(fRho, fR); }
 
//.............................................................................
void PMNS_OQS::UpdateRho() { get_SU3(fR, fRho); }
 
//.............................................................................
void PMNS_OQS::RotateState(bool to_mass)
{
  BuildHms();
  if (!to_mass) UpdateRho();
  PMNS_DensityMatrix::RotateState(to_mass, fUM);
  if (to_mass) BuildR();
}
 
//.............................................................................
void PMNS_OQS::Propagate()
{
  RotateState(true);
  PMNS_Base::Propagate();
  RotateState(false);
}
 
//.............................................................................
void PMNS_OQS::PropagatePath(NuPath p)
{
  BuildM(p);
 
  // Solve for eigensystem
  SolveHam();
 
  // Convert path length to eV
  double lengthIneV = kKm2eV * p.length;
 
  fM *= lengthIneV;
  fM = fM.exp();
 
  fRt[0] = fR[0];
 
  for (int i = 1; i < SU3_DIM; ++i) {
    fRt[i] = 0;
    for (int j = 1; j < SU3_DIM; ++j) { fRt[i] += fM(i - 1, j - 1) * fR[j]; }
  }
  fR = fRt;
}
 
//.............................................................................
matrixD PMNS_OQS::ProbMatrix(int nflvi, int nflvf)
{
  THROW_ON_INVALID_ARG(nflvi >= 0, nflvi);
  THROW_ON_INVALID_ARG(nflvi <= fNumNus, nflvi, fNumNus);
  THROW_ON_INVALID_ARG(nflvf >= 0, nflvf);
  THROW_ON_INVALID_ARG(nflvf <= fNumNus, nflvf, fNumNus);
 
  // Output probabilities
  matrixD probs(nflvi, vectorD(nflvf));
 
  // List of states
  vector<vectorD> allstates(nflvi, vectorD(9));
 
  // Reset all initial states
  for (int i = 0; i < nflvi; i++) {
    ResetToFlavour(i);
    RotateState(true);
    allstates[i] = fR;
  }
 
  // Propagate all states in parallel
  for (int i = 0; i < int(fNuPaths.size()); i++) {
    for (int flvi = 0; flvi < nflvi; flvi++) {
      fR = allstates[flvi];
      PropagatePath(fNuPaths[i]);
      allstates[flvi] = fR;
    }
  }
 
  // Get all probabilities
  for (int flvi = 0; flvi < nflvi; flvi++) {
    fR = allstates[flvi];
    RotateState(false);
    for (int flvj = 0; flvj < nflvf; flvj++) { probs[flvi][flvj] = P(flvj); }
  }
 
  return probs;
}