PMNS_LIV.cxx

Go to the documentation of this file.

 
 
#include <iostream>
 
#include "PMNS_LIV.h"
 
using namespace OscProb;
 
using namespace std;
 
//.............................................................................
PMNS_LIV::PMNS_LIV() : PMNS_Fast()
{
  SetIsOscProbAvg(true);
  SetStdPath();
 
  for (int dim = 3; dim < 8; dim += 2) {
    for (int flvi = 0; flvi < 3; flvi++) {
      for (int flvj = flvi; flvj < 3; flvj++) {
        SetaT(flvi, flvj, dim, 0, 0);
        SetcT(flvi, flvj, dim + 1, 0, 0);
      }
    }
  }
}
 
//.............................................................................
PMNS_LIV::~PMNS_LIV() {}
 
//.............................................................................
void PMNS_LIV::SetaT(int flvi, int flvj, int dim, double val, double phase)
{
  if (flvi > flvj) {
    cerr << "WARNING: First argument should be smaller or equal to second "
            "argument"
         << endl
         << "Setting reverse order (aT_" << flvj << flvi << "). " << endl;
    int temp = flvi;
    flvi     = flvj;
    flvj     = temp;
  }
  if (flvi < 0 || flvi > 2 || flvj < flvi || flvj > 2) {
    cerr << "WARNING: aT_" << flvi << flvj << " not valid for " << fNumNus
         << " neutrinos. Doing nothing." << endl;
    return;
  }
  if (dim != 3 && dim != 5 && dim != 7) {
    cerr << "WARNING: Invalid aT coefficient dimension dim=" << dim
         << " not in  [3,5,7].\n";
    return;
  }
 
  int pos = (dim - 3) / 2;
 
  complexD h = val;
 
  if (flvi != flvj) h *= complexD(cos(phase), sin(phase));
 
  bool isSame = (faT[flvi][flvj][pos] == h);
 
  if (!isSame) ClearCache();
 
  fGotES *= isSame;
 
  faT[flvi][flvj][pos] = h;
}
 
//.............................................................................
void PMNS_LIV::SetcT(int flvi, int flvj, int dim, double val, double phase)
{
  if (flvi > flvj) {
    cerr << "WARNING: First argument should be smaller or equal to second "
            "argument"
         << endl
         << "Setting reverse order (cT_" << flvj << flvi << "). " << endl;
    int temp = flvi;
    flvi     = flvj;
    flvj     = temp;
  }
  if (flvi < 0 || flvi > 2 || flvj < flvi || flvj > 2) {
    cerr << "WARNING: cT_" << flvi << flvj << " not valid for " << fNumNus
         << " neutrinos. Doing nothing." << endl;
    return;
  }
  if (dim != 4 && dim != 6 && dim != 8) {
    cerr << "WARNING: Invalid cT coefficient dimension dim=" << dim
         << " not in  [4,6,8].\n";
    return;
  }
 
  int pos = dim / 2 - 2;
 
  complexD h = val;
 
  if (flvi != flvj) h *= complexD(cos(phase), sin(phase));
 
  bool isSame = (fcT[flvi][flvj][pos] == h);
 
  if (!isSame) ClearCache();
 
  fGotES *= isSame;
 
  fcT[flvi][flvj][pos] = h;
}
 
//.............................................................................
complexD PMNS_LIV::GetaT(int flvi, int flvj, int dim)
{
  if (flvi > flvj) {
    cerr << "WARNING: First argument should be smaller or equal to second "
            "argument"
         << endl
         << "Setting reverse order (aT_" << flvj << flvi << "). " << endl;
    int temp = flvi;
    flvi     = flvj;
    flvj     = temp;
  }
  if (flvi < 0 || flvi > 2 || flvj < flvi || flvj > 2) {
    cerr << "WARNING: aT_" << flvi << flvj << " not valid for " << fNumNus
         << " neutrinos. Returning 0." << endl;
    return zero;
  }
  if (dim != 3 && dim != 5 && dim != 7) {
    cerr << "WARNING: Invalid aT coefficient dimension dim=" << dim
         << " not in  [3,5,7].\n";
    return zero;
  }
 
  int pos = (dim - 3) / 2;
 
  return faT[flvi][flvj][pos];
}
 
//.............................................................................
complex<double> PMNS_LIV::GetcT(int flvi, int flvj, int dim)
{
  if (flvi > flvj) {
    cerr << "WARNING: First argument should be smaller or equal to second "
            "argument"
         << endl
         << "Setting reverse order (cT_" << flvj << flvi << "). " << endl;
    int temp = flvi;
    flvi     = flvj;
    flvj     = temp;
  }
  if (flvi < 0 || flvi > 2 || flvj < flvi || flvj > 2) {
    cerr << "WARNING: cT_" << flvi << flvj << " not valid for " << fNumNus
         << " neutrinos. Returning 0." << endl;
    return zero;
  }
  if (dim != 4 && dim != 6 && dim != 8) {
    cerr << "WARNING: Invalid cT coefficient dimension dim=" << dim
         << " not in  [4,6,8].\n";
    return zero;
  }
 
  int pos = dim / 2 - 2;
 
  return fcT[flvi][flvj][pos];
}
 
//.............................................................................
void PMNS_LIV::SolveHam() { PMNS_Fast::SolveHamMatter(); }
 
//.............................................................................
void PMNS_LIV::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;
 
  // Add LIV terms for all dimensions
  double energy_pow = kGeV2eV;
 
  for (int dim = 3; dim < 9; dim++) {
    // Increase energy power
    if (dim > 3) energy_pow *= fEnergy;
 
    // Set parameter index
    int coeff_idx = (dim - 3) / 2;
 
    for (int i = 0; i < fNumNus; i++) {
      for (int j = i; j < fNumNus; j++) {
        complexD liv_term = energy_pow;
 
        if (dim % 2 == 1) { // aT is CPT-odd
          if (fIsNuBar)
            liv_term *= -faT[i][j][coeff_idx];
          else
            liv_term *= faT[i][j][coeff_idx];
        }
        else { // cT is CPT-even
          // Add 4/3 coefficient for backward compatibility
          if (dim == 4)
            liv_term *= -4 / 3. * fcT[i][j][coeff_idx];
          else
            liv_term *= -fcT[i][j][coeff_idx];
        }
 
        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]); }
    }
  }
}
 
//
// Obsolete functions for backward compatibility...
//
 
//.............................................................................
void PMNS_LIV::SetLIV(double aT_ee, double aT_mumu, double aT_tautau,
                      double aT_emu, double aT_etau, double aT_mutau,
                      double cT_ee, double cT_mumu, double cT_tautau,
                      double cT_emu, double cT_etau, double cT_mutau,
                      double delta_aT_emu, double delta_aT_etau,
                      double delta_aT_mutau, double delta_cT_emu,
                      double delta_cT_etau, double delta_cT_mutau)
{
  SetaT(0, 0, aT_ee, 0);
  SetaT(1, 1, aT_mumu, 0);
  SetaT(2, 2, aT_tautau, 0);
 
  SetaT(0, 1, aT_emu, delta_aT_emu);
  SetaT(0, 2, aT_etau, delta_aT_etau);
  SetaT(1, 2, aT_mutau, delta_aT_mutau);
 
  SetcT(0, 0, cT_ee, 0);
  SetcT(1, 1, cT_mumu, 0);
  SetcT(2, 2, cT_tautau, 0);
 
  SetcT(0, 1, cT_emu, delta_cT_emu);
  SetcT(0, 2, cT_etau, delta_cT_etau);
  SetcT(1, 2, cT_mutau, delta_cT_mutau);
}
 
//.............................................................................
void PMNS_LIV::SetaT(int flvi, int flvj, double val, double phase)
{
  SetaT(flvi, flvj, 3, val, phase);
}
 
//.............................................................................
void PMNS_LIV::SetcT(int flvi, int flvj, double val, double phase)
{
  SetcT(flvi, flvj, 4, val, phase);
}
 
//.............................................................................
void PMNS_LIV::SetaT_ee(double a) { SetaT(0, 0, a, 0); }
 
//.............................................................................
void PMNS_LIV::SetaT_mumu(double a) { SetaT(1, 1, a, 0); }
 
//.............................................................................
void PMNS_LIV::SetaT_tautau(double a) { SetaT(2, 2, a, 0); }
 
//.............................................................................
void PMNS_LIV::SetaT_emu(double a, double phi) { SetaT(0, 1, a, phi); }
 
//.............................................................................
void PMNS_LIV::SetaT_etau(double a, double phi) { SetaT(0, 2, a, phi); }
 
//.............................................................................
void PMNS_LIV::SetaT_mutau(double a, double phi) { SetaT(1, 2, a, phi); }
 
//.............................................................................
void PMNS_LIV::SetcT_ee(double a) { SetcT(0, 0, a, 0); }
 
//.............................................................................
void PMNS_LIV::SetcT_mumu(double a) { SetcT(1, 1, a, 0); }
 
//.............................................................................
void PMNS_LIV::SetcT_tautau(double a) { SetcT(2, 2, a, 0); }
 
//.............................................................................
void PMNS_LIV::SetcT_emu(double a, double phi) { SetcT(0, 1, a, phi); }
 
//.............................................................................
void PMNS_LIV::SetcT_etau(double a, double phi) { SetcT(0, 2, a, phi); }
 
//.............................................................................
void PMNS_LIV::SetcT_mutau(double a, double phi) { SetcT(1, 2, a, phi); }