PMNS_NSI.cxx

Go to the documentation of this file.

 
//
// Implementation of oscillations of neutrinos in matter in a
// three-neutrino framework with NSI.
//
// This  class inherits from the PMNS_Fast class
//
// jcoelho\@apc.in2p3.fr
 
#include <iostream>
 
#include "PMNS_NSI.h"
 
using namespace OscProb;
 
using namespace std;
 
//.............................................................................
PMNS_NSI::PMNS_NSI() : PMNS_Fast(), fEps()
{
  SetStdPath();
  SetNSI(0., 0., 0., 0., 0., 0., 0., 0., 0.);
  SetFermCoup(1, 0, 0);
}
 
//.............................................................................
PMNS_NSI::~PMNS_NSI() {}
 
//.............................................................................
void PMNS_NSI::SetNSI(double eps_ee, double eps_emu, double eps_etau,
                      double eps_mumu, double eps_mutau, double eps_tautau,
                      double delta_emu, double delta_etau, double delta_mutau)
{
  SetEps(0, 0, eps_ee, 0);
  SetEps(1, 1, eps_mumu, 0);
  SetEps(2, 2, eps_tautau, 0);
 
  SetEps(0, 1, eps_emu, delta_emu);
  SetEps(0, 2, eps_etau, delta_etau);
  SetEps(1, 2, eps_mutau, delta_mutau);
}
 
//.............................................................................
void PMNS_NSI::SetEps(int flvi, int flvj, double val, double phase)
{
  if (flvi > flvj) {
    cerr << "WARNING: First argument should be smaller or equal to second "
            "argument"
         << endl
         << "Setting reverse order (Eps_" << flvj << flvi << "). " << endl;
    int temp = flvi;
    flvi     = flvj;
    flvj     = temp;
  }
  if (flvi < 0 || flvi > 2 || flvj < flvi || flvj > 2) {
    cerr << "WARNING: Eps_" << flvi << flvj << " not valid for " << fNumNus
         << " neutrinos. Doing nothing." << endl;
    return;
  }
 
  complexD h = val;
 
  if (flvi != flvj) h *= complexD(cos(phase), sin(phase));
 
  bool isSame = (fEps[flvi][flvj] == h);
 
  if (!isSame) ClearCache();
 
  fGotES *= isSame;
 
  fEps[flvi][flvj] = h;
}
 
//.............................................................................
complexD PMNS_NSI::GetEps(int flvi, int flvj)
{
  if (flvi > flvj) {
    cerr << "WARNING: First argument should be smaller or equal to second "
            "argument"
         << endl
         << "Setting reverse order (Eps_" << flvj << flvi << "). " << endl;
    int temp = flvi;
    flvi     = flvj;
    flvj     = temp;
  }
  if (flvi < 0 || flvi > 2 || flvj < flvi || flvj > 2) {
    cerr << "WARNING: Eps_" << flvi << flvj << " not valid for " << fNumNus
         << " neutrinos. Returning 0." << endl;
    return zero;
  }
 
  return fEps[flvi][flvj];
}
 
//.............................................................................
void PMNS_NSI::SetEps_ee(double a) { SetEps(0, 0, a, 0); }
 
//.............................................................................
void PMNS_NSI::SetEps_mumu(double a) { SetEps(1, 1, a, 0); }
 
//.............................................................................
void PMNS_NSI::SetEps_tautau(double a) { SetEps(2, 2, a, 0); }
 
//.............................................................................
void PMNS_NSI::SetEps_emu(double a, double phi) { SetEps(0, 1, a, phi); }
 
//.............................................................................
void PMNS_NSI::SetEps_etau(double a, double phi) { SetEps(0, 2, a, phi); }
 
//.............................................................................
void PMNS_NSI::SetEps_mutau(double a, double phi) { SetEps(1, 2, a, phi); }
 
//.............................................................................
void PMNS_NSI::UpdateHam()
{
  double lv = 2 * kGeV2eV * fEnergy; // 2*E in eV
 
  double kr2GNe   = kK2 * M_SQRT2 * kGf * fPath.density;
  double kr2GNnsi = kr2GNe;
 
  kr2GNe *= fPath.zoa;      // Std matter potential in eV
  kr2GNnsi *= GetZoACoup(); // NSI matter potential in eV
 
  // Finish build Hamiltonian in matter with dimension of eV
  for (int i = 0; i < fNumNus; i++) {
    for (int j = i; j < fNumNus; j++) {
      if (!fIsNuBar)
        fHam[i][j] = fHms[i][j] / lv + kr2GNnsi * fEps[i][j];
      else
        fHam[i][j] = conj(fHms[i][j] / lv - kr2GNnsi * fEps[i][j]);
    }
  }
 
  if (!fIsNuBar)
    fHam[0][0] += kr2GNe;
  else
    fHam[0][0] -= kr2GNe;
}
 
//.............................................................................
void PMNS_NSI::SetCoupByIndex(double c, int i)
{
  bool isSame = (fNSIcoup[i] == c);
 
  if (!isSame) ClearCache();
 
  fGotES *= isSame;
 
  fNSIcoup[i] = c;
}
 
//.............................................................................
void PMNS_NSI::SetElecCoup(double e) { SetCoupByIndex(e, 0); }
 
//.............................................................................
void PMNS_NSI::SetUpCoup(double u) { SetCoupByIndex(u, 1); }
 
//.............................................................................
void PMNS_NSI::SetDownCoup(double d) { SetCoupByIndex(d, 2); }
 
//.............................................................................
void PMNS_NSI::SetFermCoup(double e, double u, double d)
{
  SetElecCoup(e); // electron coupling
  SetUpCoup(u);   // u-quark coupling
  SetDownCoup(d); // d-quark coupling
}
 
//.............................................................................
double PMNS_NSI::GetElecCoup() { return fNSIcoup[0]; }
 
//.............................................................................
double PMNS_NSI::GetUpCoup() { return fNSIcoup[1]; }
 
//.............................................................................
double PMNS_NSI::GetDownCoup() { return fNSIcoup[2]; }
 
//.............................................................................
double PMNS_NSI::GetZoACoup()
{
  return fNSIcoup[0] * fPath.zoa          // electrons: Z
         + fNSIcoup[1] * (1 + fPath.zoa)  // u-quarks:  A + Z
         + fNSIcoup[2] * (2 - fPath.zoa); // d-quarks: 2A - Z
}