PMNS_Iter.cxx

Go to the documentation of this file.

 
//
// Implementation of oscillations of neutrinos in matter in a
// three-neutrino framework.
//
// jcoelho\@apc.in2p3.fr
 
#include "PMNS_Iter.h"
 
using namespace OscProb;
 
//.............................................................................
PMNS_Iter::PMNS_Iter() : PMNS_Fast(), fPrec(1e-3) {}
 
//.............................................................................
PMNS_Iter::~PMNS_Iter() {}
 
//.............................................................................
void PMNS_Iter::SetIsNuBar(bool isNuBar)
{
  // Check if value is actually changing
  fBuiltHms *= (fIsNuBar == isNuBar);
 
  fIsNuBar = isNuBar;
}
 
//.............................................................................
void PMNS_Iter::SetExpVL(NuPath p)
{
  double kr2GNe = kK2 * M_SQRT2 * kGf;
  kr2GNe *= p.density * p.zoa; // Matter potential in eV
 
  fVL = kr2GNe * kKm2eV * p.length;
 
  fExpVL = complexD(cos(fVL), -sin(fVL));
  if (fIsNuBar) fExpVL = conj(fExpVL);
}
 
//.............................................................................
void PMNS_Iter::PropMatter() { fNuState[0] *= fExpVL; }
 
//.............................................................................
void PMNS_Iter::SolveHam()
{
  // Do vacuum oscillation
  if (!fBuiltHms) {
    PMNS_Fast::SetVacuumEigensystem();
    fBuiltHms   = true;
    fGotES      = true;
    fPrevEnergy = fEnergy;
    return;
  }
 
  if (fGotES) return;
 
  for (int i = 1; i < fNumNus; i++) { fEval[i] *= fPrevEnergy / fEnergy; }
 
  fPrevEnergy = fEnergy;
 
  fGotES = true;
  return;
}
 
//.............................................................................
void PMNS_Iter::SetPrec(double prec)
{
  fPrec = prec;
  if (fPrec <= 0) fPrec = 1e-3;
}
 
//.............................................................................
void PMNS_Iter::PropagatePath(NuPath p)
{
  SetExpVL(p);
 
  double dm = 0;
  for (int i = 0; i < fNumNus; i++) {
    if (dm < fabs(fDm[i])) dm = fabs(fDm[i]);
  }
  dm *= kKm2eV * p.length / (2 * kGeV2eV * fEnergy);
 
  int nsplit = sqrt(0.065 * dm * fVL / fPrec);
  if (nsplit > 1) {
    p.length /= nsplit;
    SetExpVL(p);
  }
  else
    nsplit = 1;
 
  for (int i = 0; i < nsplit; i++) {
    PMNS_Base::PropagatePath(p);
    PropMatter();
  }
}