/////////////////////////////////////////////////////////////////////////////// // // Implementation of oscillations of neutrinos in matter in a // three-neutrino framework. // // jcoelho\@apc.in2p3.fr /////////////////////////////////////////////////////////////////////////////// #include "PMNS_Iter.h" using namespace OscProb; //............................................................................. /// /// Constructor. \sa PMNS_Fast::PMNS_Fast /// /// This class is restricted to 3 neutrino flavours. /// PMNS_Iter::PMNS_Iter() : PMNS_Fast(), fPrec(1e-3) {} //............................................................................. /// /// Nothing to clean. /// PMNS_Iter::~PMNS_Iter() {} //............................................................................. /// /// Propagate neutrino through matter component. /// 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); } //............................................................................. /// /// Propagate neutrino through matter component. /// void PMNS_Iter::PropMatter() { fNuState[0] *= fExpVL; } //............................................................................. /// /// Just use the vacuum Hamiltonian to start. /// void PMNS_Iter::SolveHam() { if (fGotES) return; // Do vacuum oscillation if (!fBuiltHms) { fPath.density = 0; PMNS_Fast::SolveHam(); fBuiltHms = true; } else { for (int i = 1; i < fNumNus; i++) { fEval[i] *= fPrevEnergy / fEnergy; } } fPrevEnergy = fEnergy; fGotES = true; return; } //............................................................................. /// /// Set iterator precision /// void PMNS_Iter::SetPrec(double prec) { fPrec = prec; if (fPrec <= 0) fPrec = 1e-3; } //............................................................................. /// /// Propagate neutrino state through split path /// 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(); } } ///////////////////////////////////////////////////////////////////////////////