PMNS_SiderealLIV.cxx

Go to the documentation of this file.

 
 
#include <algorithm>
#include <cmath>
#include <iostream>
 
#include "MatrixDecomp/zheevh3.h"
#include "PMNS_SiderealLIV.h"
 
using namespace OscProb;
 
using namespace std;
 
//.............................................................................
PMNS_SiderealLIV::PMNS_SiderealLIV()
    : PMNS_Fast(), fa{}, fc{}, fN{0.0, 0.0, 0.0}, fTime(0.0), fChi(0.0)
{
  SetStdPath();
  // The Maltoni averaging method assumes H ∝ 1/E, which is not valid here
  // because the LIV terms scale as E^0 (aT) and E^1 (cT).
  SetIsOscProbAvg(true);
}
 
//.............................................................................
PMNS_SiderealLIV::~PMNS_SiderealLIV() {}
 
//.............................................................................
bool PMNS_SiderealLIV::ValidateAndOrder(int& flvi, int& flvj, const char* name)
{
  if (flvi > flvj) {
    cerr << "WARNING: First argument should be smaller or equal to second "
            "argument"
         << endl
         << "Setting reverse order (" << name << "_" << flvj << flvi << "). "
         << endl;
    std::swap(flvi, flvj);
  }
  if (flvi < 0 || flvi > 2 || flvj < flvi || flvj > 2) {
    cerr << "WARNING: " << name << "_" << flvi << flvj
         << " not valid for 3 neutrinos. Doing nothing." << endl;
    return false;
  }
  return true;
}
 
//.............................................................................
bool PMNS_SiderealLIV::ValidateCoord(int coord, const char* coordName,
                                     const char* paramName)
{
  if (coord < 0 || coord > 2) {
    cerr << "WARNING: " << coordName << " = " << coord
         << " out of range [0,2] for " << paramName << ". Doing nothing."
         << endl;
    return false;
  }
  return true;
}
 
//.............................................................................
void PMNS_SiderealLIV::SetA(int flvi, int flvj, int coord, double val)
{
  if (!ValidateAndOrder(flvi, flvj, "aT")) return;
  if (!ValidateCoord(coord, "coord", "aT")) return;
 
  fa[flvi][flvj][coord] = val;
}
 
//.............................................................................
void PMNS_SiderealLIV::SetC(int flvi, int flvj, int coord1, int coord2,
                            double val)
{
  if (!ValidateAndOrder(flvi, flvj, "cT")) return;
  if (!ValidateCoord(coord1, "coord1", "cT")) return;
  if (!ValidateCoord(coord2, "coord2", "cT")) return;
 
  fc[flvi][flvj][coord1][coord2] = val;
}
 
//.............................................................................
double PMNS_SiderealLIV::GetA(int flvi, int flvj, int coord)
{
  if (!ValidateAndOrder(flvi, flvj, "aT")) return 0;
  if (!ValidateCoord(coord, "coord", "aT")) return 0;
 
  return fa[flvi][flvj][coord];
}
 
//.............................................................................
double PMNS_SiderealLIV::GetC(int flvi, int flvj, int coord1, int coord2)
{
  if (!ValidateAndOrder(flvi, flvj, "cT")) return 0;
  if (!ValidateCoord(coord1, "coord1", "cT")) return 0;
  if (!ValidateCoord(coord2, "coord2", "cT")) return 0;
 
  return fc[flvi][flvj][coord1][coord2];
}
 
//.............................................................................
void PMNS_SiderealLIV::SetNeutrinoDirection(double zenith, double azimuth)
{
  double chi = fChi * M_PI / 180.0;
  double zen = zenith * M_PI / 180.0;
  double azi = azimuth * M_PI / 180.0;
  double N0  = cos(chi) * sin(zen) * cos(azi) + sin(chi) * cos(zen);
  double N1  = sin(zen) * sin(azi);
  double N2  = -sin(chi) * sin(zen) * cos(azi) + cos(chi) * cos(zen);
 
  fGotES *= (fN[0] == N0 && fN[1] == N1 && fN[2] == N2);
 
  fN[0] = N0;
  fN[1] = N1;
  fN[2] = N2;
}
 
//.............................................................................
void PMNS_SiderealLIV::SetTimeHours(double hours)
{
  fGotES *= (fTime == hours);
  fTime = hours;
}
 
//.............................................................................
void PMNS_SiderealLIV::SetColatitude(double chi) { fChi = chi; }
 
//.............................................................................
void PMNS_SiderealLIV::SetColatitude(double deg, double min, double sec)
{
  double latitude = deg + min / 60.0 + sec / 3600.0;
  fChi            = 90.0 - latitude;
}
 
//.............................................................................
double PMNS_SiderealLIV::GetColatitude() const { return fChi; }
 
//.............................................................................
void PMNS_SiderealLIV::SolveHam() { PMNS_Fast::SolveHamMatter(); }
 
//.............................................................................
void PMNS_SiderealLIV::UpdateHam()
{
  double lv = 2 * kGeV2eV * fEnergy; // 2*E in eV
 
  // Set the vacuum Hamiltonian
  for (int i = 0; i < fNumNus; i++) {
    for (int j = i; j < fNumNus; j++) { fHam[i][j] = fHms[i][j] / lv; }
  }
 
  // Add matter potential
  double kr2GNe = kK2 * M_SQRT2 * kGf * fPath.density * fPath.zoa;
  if (!fIsNuBar)
    fHam[0][0] += kr2GNe;
  else
    fHam[0][0] -= kr2GNe;
 
  for (int i = 0; i < fNumNus; i++) {
    for (int j = i; j < fNumNus; j++) {
      double sign = fIsNuBar ? -kGeV2eV : kGeV2eV;
 
      double C0 = -sign * fa[i][j][2] * fN[2];
 
      double As0 = sign * (fa[i][j][0] * fN[1] - fa[i][j][1] * fN[0]);
      double Ac0 = -sign * (fa[i][j][0] * fN[0] + fa[i][j][1] * fN[1]);
 
      double As1 = kGeV2eV * (2 * fN[1] * fN[2] * fc[i][j][0][2] -
                               2 * fN[0] * fN[2] * fc[i][j][1][2]);
      double Ac1 = kGeV2eV * (-2 * fN[0] * fN[2] * fc[i][j][0][2] -
                               2 * fN[1] * fN[2] * fc[i][j][1][2]);
 
      double Bs1 =
          kGeV2eV * (fN[0] * fN[1] * (fc[i][j][0][0] - fc[i][j][1][1]) -
                     (fN[0] * fN[0] - fN[1] * fN[1]) * fc[i][j][0][1]);
      double Bs = fEnergy * Bs1;
 
      double Bc1 = kGeV2eV * (-0.5 * (fN[0] * fN[0] - fN[1] * fN[1]) *
                                   (fc[i][j][0][0] - fc[i][j][1][1]) -
                               2.0 * fN[0] * fN[1] * fc[i][j][0][1]);
      double Bc = fEnergy * Bc1;
 
      double As = As0 + fEnergy * As1;
      double Ac = Ac0 + fEnergy * Ac1;
 
      double liv_term =
          C0 + As * sin(kOmegaSidereal * fTime) +
          Ac * cos(kOmegaSidereal * fTime) +
          Bs * 2 * sin(kOmegaSidereal * fTime) * cos(kOmegaSidereal * fTime) +
          Bc * (cos(kOmegaSidereal * fTime) * cos(kOmegaSidereal * fTime) -
                sin(kOmegaSidereal * fTime) * sin(kOmegaSidereal * fTime));
 
      fHam[i][j] += liv_term;
    }
  }
 
  // Conjugate Hamiltonian for antineutrinos
  if (fIsNuBar) {
    for (int i = 0; i < fNumNus; i++) {
      for (int j = i + 1; j < fNumNus; j++) { fHam[i][j] = conj(fHam[i][j]); }
    }
  }
}