Oscillation probability calculators. More...

Classes
class	_IOscCalculator
	General interface to oscillation calculators. More...

class	_IOscCalculatorAdjustable
	General interface to any calculator that lets you set the parameters. More...

class	_NoOscillations
	Pass neutrinos through unchanged. More...

struct	_OscCache

class	_OscCalculatorAnalytic

class	_OscCalculatorDMP
	Helper struct for the cache. Might not need this. More...

class	_OscCalculatorPMNS
	Adapt the PMNS calculator to standard interface. More...

class	_OscCalculatorPMNSOpt
	Optimized version of OscCalculatorPMNS. More...

struct	_OscParameters

class	_PMNS

class	_PMNS_DMP

class	_PMNSOpt
	Optimized version of PMNS. More...

struct	cmplx
	std::complex takes a lot of care with inf/nan which we don't want More...

class	EarthModel

struct	EigenSystem
	Helper struct for the cache. Might not need this. More...

class	OscCalculator

class	OscCalculatorCPT

class	OscCalculatorDumb
	Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it's some kind of average of all of those. More...

class	OscCalculatorGeneral
	More generic (but probably slower) oscillation calculations. More...

class	OscCalculatorPMNS_CPT
	Adapt the PMNS calculator to standard interface and include neutrino and anti neutrino oscillations seperately. More...

class	OscCalculatorPMNS_NSI
	Optimized version of OscCalculatorPMNS. More...

class	OscCalculatorPMNSOptEigen
	A re-optimized version of OscCalculatorPMNSOpt. More...

class	OscCalculatorSterile
	Adapt the PMNS_Sterile calculator to standard interface. More...

class	OscCalculatorSterileBeam
	Adapt the PMNS_Sterile calculator to standard interface. More...

class	OscCalculatorSterileTrivial
	Version of OscCalculatorSterile that always returns probability of 1. More...

class	PMNS_NSI

class	PMNS_Sterile

Typedefs
typedef _IOscCalculator< double >	IOscCalculator

typedef _IOscCalculatorAdjustable< double >	IOscCalculatorAdjustable

typedef _IOscCalculatorAdjustable< stan::math::var >	IOscCalculatorAdjustableStan

typedef _IOscCalculator< stan::math::var >	IOscCalculatorStan

typedef _OscCalculatorDMP< stan::math::var >	OscCalculatorDMPStan

typedef _OscCalculatorPMNS< stan::math::var >	OscCalculatorPMNSStan

typedef _OscCalculatorPMNSOpt< stan::math::var >	OscCalculatorPMNSOptStan

typedef _NoOscillations< double >	NoOscillations

typedef _OscCalculatorAnalytic< double >	OscCalculatorAnalytic

using	SDMap = std::map< const ana::SigmaDelta *, std::pair< double, double > >

typedef _OscCalculatorDMP< double >	OscCalculatorDMP

typedef std::complex< long double >	val_t

typedef ublas::bounded_array< val_t, kNumFlavours >	alloc_t

typedef ublas::bounded_matrix< val_t, kNumFlavours, kNumFlavours >	ComplexMat

typedef ublas::c_vector< val_t, kNumFlavours >	ComplexVec

typedef ublas::unit_vector< val_t, alloc_t >	UnitVec

typedef _OscCalculatorPMNS< double >	OscCalculatorPMNS

typedef _OscCalculatorPMNSOpt< double >	OscCalculatorPMNSOpt

typedef _OscParameters< double >	OscParameters

typedef _PMNS< double >	PMNS

typedef _PMNSOpt< double >	PMNSOpt

Enumerations
enum	ENuSign { ENuSign::kNu, ENuSign::kNuBar }
	Simple neutrino/anti-neutrino enum. More...

Functions
template<typename T , typename U >
void	CopyParams (const osc::_IOscCalculatorAdjustable< T > inCalc, osc::_IOscCalculatorAdjustable< U > outCalc)
	Copy parameters from one calculator to another, irrespective of their type. More...

template void	CopyParams (const osc::_IOscCalculatorAdjustable< double > inCalc, osc::_IOscCalculatorAdjustable< stan::math::var > outCalc)

template void	CopyParams (const osc::_IOscCalculatorAdjustable< stan::math::var > inCalc, osc::_IOscCalculatorAdjustable< double > outCalc)

template void	CopyParams (const osc::_IOscCalculatorAdjustable< stan::math::var > inCalc, osc::_IOscCalculatorAdjustable< stan::math::var > outCalc)

template void	CopyParams (const osc::_IOscCalculatorAdjustable< double > inCalc, osc::_IOscCalculatorAdjustable< double > outCalc)

template<class T >
	__attribute__ ((always_inline)) cmplx< T > operator*(const cmplx< T > &x

template<class T , class U >
	__attribute__ ((always_inline)) cmplx< T > operator*(const U &x

template<class T >
T	sqr (T x)

template<class T >
T	cube (T x)

template<class T >
std::array< T, 3 >	SolveCubic (T b, T c, T d)
	Solve x^3 + bx^2 + cx + d = 0. More...

const OscCalculatorCPT *	DowncastToCPT (const IOscCalculatorAdjustable *osc)

OscCalculatorCPT *	DowncastToCPT (IOscCalculatorAdjustable *osc)

void	_sincos (double theta, double &s, double &c)

const ublas::zero_matrix< val_t, alloc_t >	kZeroMat (kNumFlavours, kNumFlavours)

const ublas::identity_matrix< val_t, alloc_t >	kIdentity (kNumFlavours)

ComplexMat	GetPMNS (OscCalculatorGeneral::Priv *d)

ComplexMat	VacuumHamiltonian (const ComplexMat &U, const std::vector< long double > &mSq, long double E)

ComplexMat	MatterHamiltonianComponent (long double Ne, long double emutau)

ComplexMat	MatrixExp (const ComplexMat &m2)

ComplexVec	EvolveState (ComplexVec A, const ComplexMat &H, long double L)

void	conjugate_elements (ComplexMat &m)

const OscCalculatorPMNS_NSI *	DowncastToNSI (const IOscCalculator *calc)

OscCalculatorPMNS_NSI *	DowncastToNSI (IOscCalculator *calc)

Eigen::Array33d	Probabilities (const Eigen::Matrix3cd &to, const Eigen::Matrix3cd &from)

const OscCalculatorSterile *	DowncastToSterile (const IOscCalculator *calc)

OscCalculatorSterile *	DowncastToSterile (IOscCalculator *calc)

const OscCalculatorSterileBeam *	DowncastToSterileBeam (const IOscCalculator *calc)

OscCalculatorSterileBeam *	DowncastToSterileBeam (IOscCalculator *calc)

template<typename T >
std::ostream &	operator<< (std::ostream &os, _OscParameters< T > const &p)

template<typename T >
std::ostream &	operator<< (std::ostream &os, _PMNS_DMP< T > &p)

Variables
const cmplx< T > &	y

const double	sqrt3 = sqrt(3)

const unsigned int	kNumFlavours = 3

static double	eVsqkm_to_GeV = 1e-9 / 1.973269681602260e-7 * 1e3

static double	YerhoE2a = 1.52e-4

static const double	kKm2eV = 5.06773103202e+09
	km to eV^-1 More...

static const double	kK2 = 4.62711492217e-09
	mole/GeV^2/cm^3 to eV More...

static const double	kGeV2eV = 1.0e+09
	GeV to eV. More...

static const double	kGf = 1.166371e-5

Detailed Description

Oscillation probability calculators.

Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it's some kind of average of all of those.

Typedef Documentation

typedef ublas::bounded_array<val_t, kNumFlavours> osc::alloc_t

Definition at line 27 of file OscCalculatorGeneral.cxx.

typedef ublas::bounded_matrix<val_t, kNumFlavours, kNumFlavours> osc::ComplexMat

Definition at line 29 of file OscCalculatorGeneral.cxx.

typedef ublas::c_vector<val_t, kNumFlavours> osc::ComplexVec

Definition at line 30 of file OscCalculatorGeneral.cxx.

typedef _IOscCalculator< double > osc::IOscCalculator

Definition at line 16 of file FwdDeclare.h.

typedef _IOscCalculatorAdjustable< double > osc::IOscCalculatorAdjustable

Definition at line 19 of file FwdDeclare.h.

typedef _IOscCalculatorAdjustable<stan::math::var> osc::IOscCalculatorAdjustableStan

Definition at line 19 of file StanTypedefs.h.

typedef _IOscCalculator<stan::math::var> osc::IOscCalculatorStan

Definition at line 22 of file StanTypedefs.h.

typedef _NoOscillations<double> osc::NoOscillations

Definition at line 62 of file IOscCalculator.h.

typedef _OscCalculatorAnalytic<double> osc::OscCalculatorAnalytic

Definition at line 104 of file OscCalculatorAnalytic.h.

typedef _OscCalculatorDMP<double> osc::OscCalculatorDMP

Definition at line 81 of file OscCalculatorDMP.h.

typedef _OscCalculatorDMP<stan::math::var> osc::OscCalculatorDMPStan

Definition at line 25 of file StanTypedefs.h.

typedef _OscCalculatorPMNS<double> osc::OscCalculatorPMNS

Definition at line 57 of file OscCalculatorPMNS.h.

typedef _OscCalculatorPMNSOpt<double> osc::OscCalculatorPMNSOpt

Definition at line 72 of file OscCalculatorPMNSOpt.h.

typedef _OscCalculatorPMNSOpt<stan::math::var> osc::OscCalculatorPMNSOptStan

Definition at line 31 of file StanTypedefs.h.

typedef _OscCalculatorPMNS<stan::math::var> osc::OscCalculatorPMNSStan

Definition at line 28 of file StanTypedefs.h.

typedef _OscParameters<double> osc::OscParameters

Definition at line 33 of file OscParameters.h.

typedef _PMNS<double> osc::PMNS

Definition at line 226 of file PMNS.h.

typedef _PMNSOpt<double> osc::PMNSOpt

Definition at line 130 of file PMNSOpt.h.

using osc::SDMap = typedef std::map< const ana::SigmaDelta*, std::pair<double, double> >

Definition at line 33 of file OscCalculatorCPT.h.

typedef ublas::unit_vector<val_t, alloc_t> osc::UnitVec

Definition at line 31 of file OscCalculatorGeneral.cxx.

typedef std::complex<long double> osc::val_t

Definition at line 26 of file OscCalculatorGeneral.cxx.

Enumeration Type Documentation

enum osc::ENuSign

strong

Simple neutrino/anti-neutrino enum.

Enumerator
kNu
kNuBar

Definition at line 28 of file OscCalculatorCPT.h.

28 { kNu, kNuBar };

osc::ENuSign::kNuBar

ana::Sign::kNu

Neutrinos-only.

Definition: IPrediction.h:49

Function Documentation

template<class T >

osc::__attribute__ ( (always_inline) ) const

inline

Definition at line 74 of file OscCalculatorAnalytic.cxx.

Referenced by LZ4_isLittleEndian(), mouseover(), and sincos().

   {
     return cmplx<T>(-x.re, -x.im);
   }

template<class T , class U >

osc::__attribute__ ( (always_inline) ) const

inline

void osc::_sincos	(	double	theta,
		double &	s,
		double &	c
	)

inline

Definition at line 26 of file OscCalculatorDMP.cxx.

References sincos().

Referenced by osc::OscCalculatorPMNSOptEigen::BuildHam(), and osc::OscCalculatorPMNSOptEigen::SaveLastParams().

   {
 #ifdef __APPLE_CC__
     __sincos(theta, &s, &c);
 #else
     sincos(theta, &s, &c);
 #endif
   }

void osc::conjugate_elements ( ComplexMat & m )

Definition at line 227 of file OscCalculatorGeneral.cxx.

References MECModelEnuComparisons::i, calib::j, kNumFlavours, and m.

Referenced by osc::OscCalculatorGeneral::P().

   {
     for(unsigned int i = 0; i < kNumFlavours; ++i)
       for(unsigned int j = 0; j < kNumFlavours; ++j)
         m(i, j) = std::conj(m(i, j));
   }

template<typename T , typename U >

void osc::CopyParams	(	const osc::_IOscCalculatorAdjustable< T > *	inCalc,
		osc::_IOscCalculatorAdjustable< U > *	outCalc
	)

Copy parameters from one calculator to another, irrespective of their type.

Definition at line 47 of file IOscCalculator.cxx.

Referenced by ana::StanFitter::CreateCalculator(), ana::StanFitter::log_prob(), and test_stanfit_systpulls().

   {
     assert (inCalc);
     assert (outCalc);
 
     outCalc->SetL(inCalc->GetL());
     outCalc->SetRho(inCalc->GetRho());
 
     outCalc->SetdCP(util::GetValAs<U>(inCalc->GetdCP()));
     outCalc->SetDmsq21(util::GetValAs<U>(inCalc->GetDmsq21()));
     outCalc->SetDmsq32(util::GetValAs<U>(inCalc->GetDmsq32()));
     outCalc->SetTh12(util::GetValAs<U>(inCalc->GetTh12()));
     outCalc->SetTh13(util::GetValAs<U>(inCalc->GetTh13()));
     outCalc->SetTh23(util::GetValAs<U>(inCalc->GetTh23()));
   }

template void osc::CopyParams	(	const osc::_IOscCalculatorAdjustable< double > *	inCalc,
		osc::_IOscCalculatorAdjustable< stan::math::var > *	outCalc
	)

template void osc::CopyParams	(	const osc::_IOscCalculatorAdjustable< stan::math::var > *	inCalc,
		osc::_IOscCalculatorAdjustable< double > *	outCalc
	)

template void osc::CopyParams	(	const osc::_IOscCalculatorAdjustable< stan::math::var > *	inCalc,
		osc::_IOscCalculatorAdjustable< stan::math::var > *	outCalc
	)

template void osc::CopyParams	(	const osc::_IOscCalculatorAdjustable< double > *	inCalc,
		osc::_IOscCalculatorAdjustable< double > *	outCalc
	)

template<class T >

T osc::cube ( T x )

Definition at line 194 of file OscCalculatorAnalytic.cxx.

References submit_syst::x.

Referenced by SolveCubic().

194 {return x*x*x;}

submit_syst.x

x

Definition: submit_syst.py:249

const OscCalculatorCPT * osc::DowncastToCPT ( const IOscCalculatorAdjustable * osc )

Definition at line 232 of file OscCalculatorCPT.cxx.

References ana::assert().

   {
     const OscCalculatorCPT* cpt = dynamic_cast<const OscCalculatorCPT*>(osc);
     assert( cpt && "Must use OscCalculatorCPT with CPT FitVars." );
     return cpt;
   }

OscCalculatorCPT * osc::DowncastToCPT ( IOscCalculatorAdjustable * osc )

Definition at line 239 of file OscCalculatorCPT.cxx.

References ana::assert().

   {
     OscCalculatorCPT* cpt = dynamic_cast<OscCalculatorCPT*>(osc);
     assert( cpt && "Must use OscCalculatorCPT with CPT FitVars." );
     return cpt;
   }

const OscCalculatorPMNS_NSI * osc::DowncastToNSI ( const IOscCalculator * calc )

Definition at line 117 of file OscCalculatorPMNS_NSI.cxx.

References calc, om::cout, and allTimeWatchdog::endl.

Referenced by ana::FitEpsMuTau::GetValue(), and ana::FitEpsMuTau::SetValue().

   {
     const OscCalculatorPMNS_NSI* calc_nsi = dynamic_cast<const OscCalculatorPMNS_NSI*>(calc);
     if(calc_nsi) return calc_nsi;
     else             std::cout << "Input calculator was not of type OscCalculatorPMNS_NSI." << std::endl;
     return nullptr; // If the cast failed, calc_nsi should be nullptr anyway (?)                    
   }

OscCalculatorPMNS_NSI * osc::DowncastToNSI ( IOscCalculator * calc )

Definition at line 126 of file OscCalculatorPMNS_NSI.cxx.

References calc, om::cout, and allTimeWatchdog::endl.

   {
     OscCalculatorPMNS_NSI* calc_nsi = dynamic_cast<OscCalculatorPMNS_NSI*>(calc);
     if(calc_nsi) return calc_nsi;
     else             std::cout << "Input calculator was not of type OscCalculatorPMNS_NSI." << std::endl;
     return nullptr; // If the cast failed, calc_nsi should be nullptr anyway (?)                    
   }

const OscCalculatorSterile * osc::DowncastToSterile ( const IOscCalculator * calc )

Definition at line 228 of file OscCalculatorSterile.cxx.

References calc, om::cout, and allTimeWatchdog::endl.

   {
     const OscCalculatorSterileTrivial* calc_trivial
             = dynamic_cast<const OscCalculatorSterileTrivial*>(calc);
     if (calc_trivial) return calc_trivial; 
     const OscCalculatorSterile* calc_sterile = dynamic_cast<const OscCalculatorSterile*>(calc);
     if(calc_sterile) return calc_sterile;
     else             std::cout << "Input calculator was not of type OscCalculatorSterile." << std::endl;
     return nullptr; // If the cast failed, calc_sterile should be nullptr anyway
   }

OscCalculatorSterile * osc::DowncastToSterile ( IOscCalculator * calc )

Definition at line 240 of file OscCalculatorSterile.cxx.

References calc, om::cout, and allTimeWatchdog::endl.

   {
     OscCalculatorSterileTrivial* calc_trivial
             = dynamic_cast<OscCalculatorSterileTrivial*>(calc);
     if (calc_trivial) return calc_trivial;
     OscCalculatorSterile* calc_sterile = dynamic_cast<OscCalculatorSterile*>(calc);
     if(calc_sterile) return calc_sterile;
     else             std::cout << "Input calculator was not of type OscCalculatorSterile." << std::endl;
     return nullptr; // If the cast failed, calc_sterile should be nullptr anyway
   }

const OscCalculatorSterileBeam * osc::DowncastToSterileBeam ( const IOscCalculator * calc )

Definition at line 143 of file OscCalculatorSterileBeam.cxx.

References calc, om::cout, and allTimeWatchdog::endl.

   {
     //const IOscCalculator* base = new osc::OscCalculatorSterileBeam();
     const OscCalculatorSterileBeam* calc_sterile = dynamic_cast<const OscCalculatorSterileBeam*>(calc);
     if(calc_sterile) return calc_sterile;
     else             std::cout << "Input calculator was not of type OscCalculatorSterileBeam." << std::endl;
     return nullptr; // If the cast failed, calc_sterile should be nullptr anyway
   }

OscCalculatorSterileBeam * osc::DowncastToSterileBeam ( IOscCalculator * calc )

std::cout << "calc is a osc::NoOscillations calculator!" << std::endl;

Definition at line 153 of file OscCalculatorSterileBeam.cxx.

References calc, om::cout, and allTimeWatchdog::endl.

   {
     //IOscCalculator* base = new osc::OscCalculatorSterileBeam();
     //std::cout << "calc's type is: " << typeid(calc).name() << std::endl;
     //if (typeid(calc) == typeid(osc::NoOscillations))
     //  {
     /// std::cout << "calc is a osc::NoOscillations calculator!" << std::endl;
     //}
     NoOscillations* calc_noosc = dynamic_cast<NoOscillations*>(calc);
     if (calc_noosc) 
       {
         std::cout << "I was successfully cast to NoOscillations" << std::endl;
       }
 
     OscCalculatorSterileBeam* calc_sterile = dynamic_cast<OscCalculatorSterileBeam*>(calc);
     if(calc_sterile) return calc_sterile;
     else             std::cout << "Input calculator was not of type OscCalculatorSterileBeam." << std::endl;
     return nullptr; // If the cast failed, calc_sterile should be nullptr anyway
   }

ComplexVec osc::EvolveState	(	ComplexVec	A,
		const ComplexMat &	H,
		long double	L
	)

Definition at line 220 of file OscCalculatorGeneral.cxx.

References MECModelEnuComparisons::i, MatrixExp(), and stan::math::prod().

Referenced by osc::OscCalculatorGeneral::P().

   {
     const std::complex<long double> i = std::complex<long double>(0, 1);
 
     return ublas::prod(MatrixExp(-H*L*i), A);
   }

ComplexMat osc::GetPMNS ( OscCalculatorGeneral::Priv * d )

Definition at line 138 of file OscCalculatorGeneral.cxx.

References osc::OscCalculatorGeneral::Priv::atmos, osc::OscCalculatorGeneral::Priv::dirty, osc::OscCalculatorGeneral::Priv::pmns, stan::math::prod(), osc::OscCalculatorGeneral::Priv::react, and osc::OscCalculatorGeneral::Priv::solar.

Referenced by osc::OscCalculatorGeneral::P().

   {
     if(d->dirty){
       ublas::noalias(d->pmns) = ublas::prod(d->atmos, ublas::prod<ComplexMat>(d->react, d->solar));
       d->dirty = false;
     }
     return d->pmns;
   }

const ublas::identity_matrix<val_t, alloc_t> osc::kIdentity ( kNumFlavours )

Referenced by MatrixExp().

const ublas::zero_matrix<val_t, alloc_t> osc::kZeroMat	(	kNumFlavours	,
		kNumFlavours
	)

Referenced by MatterHamiltonianComponent(), and VacuumHamiltonian().

ComplexMat osc::MatrixExp ( const ComplexMat & m2 )

Definition at line 204 of file OscCalculatorGeneral.cxx.

References kIdentity(), m, getGoodRuns4SAM::n, stan::math::prod(), and runNovaSAM::ret.

Referenced by EvolveState().

   {
     // We use the identity e^(a*b) = (e^a)^b
     // First: ensure the exponent is really small, 65536 = 2^16
     ComplexMat m = m2/65536;
 
     // To first order e^x = 1+x
     ComplexMat ret = kIdentity+m;
 
     // Raise the result to the power 65536, by squaring it 16 times
     for(int n = 0; n < 16; ++n) ret = ublas::prod(ret, ret);
 
     return ret;
   }

ComplexMat osc::MatterHamiltonianComponent	(	long double	Ne,
		long double	emutau
	)

Definition at line 178 of file OscCalculatorGeneral.cxx.

References kZeroMat(), and std::sqrt().

Referenced by osc::OscCalculatorGeneral::P().

   {
     // Need to convert avogadro's constant so that the total term comes out in
     // units of inverse distance. Note that Ne will be specified in g/cm^-3
     // I put the following into Wolfram Alpha:
     // (fermi coupling constant)*((avogadro number)/cm^3)*(reduced planck constant)^2*(speed of light)^2
     // And then multiplied by 1000 because we specify L in km not m.
     const long double GF = 1.368e-4;
 
     ComplexMat H = kZeroMat;
 
     H(0, 0) = sqrt(2)*GF*Ne;
 
     // Ignoring conjugates here because we assume e_mutau is real
     H(1, 2) = H(2, 1) = emutau*sqrt(2)*GF*Ne;
 
     return H;
   }

template<typename T >

std::ostream& osc::operator<<	(	std::ostream &	os,
		_OscParameters< T > const &	p
	)

inline

Definition at line 36 of file OscParameters.h.

   {
     os << "\tdelta m^2_{21}: " << p.dmsq21
        << "\tdelta m^2_{32}: " << p.dmsq32
        << "\tdelta m^2_{31}: " << p.Dmsq31()
        << "\ttheta_{12}: " << p.th12
        << "\ttheta_{13}: " << p.th13
        << "\ttheta_{23}: " << p.th23
        << "\tdelta CP: " << p.deltacp
        << "\n";
     return os;
   }

template<typename T >

std::ostream& osc::operator<<	(	std::ostream &	os,
		_PMNS_DMP< T > &	p
	)

inline

Definition at line 103 of file PMNS_DMP.h.

     {
         os << "\tDmsqee: " << p.Dmsqee()
            //<< "\tdelta m^2_{32}: " << p.dmsq32
            //<< "\tdelta m^2_{31}: " << p.Dmsq31()
            //<< "\ttheta_{12}: " << p.th12
            //<< "\ttheta_{13}: " << p.th13
            //<< "\ttheta_{23}: " << p.th23
            //<< "\tdelta CP: " << p.deltacp
            << "\n";
       return os;
     }

Eigen::Array33d osc::Probabilities	(	const Eigen::Matrix3cd &	to,
		const Eigen::Matrix3cd &	from
	)

inline

Definition at line 28 of file OscCalculatorPMNSOptEigen.cxx.

Referenced by osc::OscCalculatorPMNSOptEigen::FillCache(), and osc::OscCalculatorPMNSOptEigen::P().

   {
     return (to * from.adjoint()).array().abs2();
   }

template<class T >

std::array<T, 3> osc::SolveCubic	(	T	b,
		T	c,
		T	d
	)

Solve x^3 + b*x^2 + c*x + d = 0.

Definition at line 200 of file OscCalculatorAnalytic.cxx.

   {
     b /= 3;
     // Now solving x^3 + 3*b*x^2 + c*x + d = 0
 
     // https://en.wikipedia.org/wiki/Cubic_equation#Depressed_cubic
     const T p = c/3 - sqr(b);
     const T q = 2*cube(b) - b*c + d;
     // Now solving t^3 + 3*p*t + q = 0
 
     // https://en.wikipedia.org/wiki/Cubic_equation#Trigonometric_solution_for_three_real_roots
     //
     // Cardano's formula looks simpler, but you still have to involve trig
     // operations anyway for the cube root, and it's harder to keep everything
     // explicitly real.
 
     const T s = sqrt(-p);
 
     const T r = acos(q/(2*p*s)) / 3;
 
     // t_k = 2*s * cos(r/3 + 2pi*k/3)
 
     // Use cos(a+b) = cosa*cosb - sina*sinb to save one trig operation
     T sinr, cosr;
     sincos(r, &sinr, &cosr);
 
     const T t0 = 2*s*cosr;
     const T t1 = s*(sqrt3*sinr - cosr);
     const T t2 = -t0-t1;
 
     return {t0 - b, t1 - b, t2 - b};
   }

template<class T >

T osc::sqr ( T x )

Definition at line 193 of file OscCalculatorAnalytic.cxx.

References submit_syst::x.

Referenced by SolveCubic().

193 {return x*x;}

submit_syst.x

x

Definition: submit_syst.py:249

ComplexMat osc::VacuumHamiltonian	(	const ComplexMat &	U,
		const std::vector< long double > &	mSq,
		long double	E
	)

Definition at line 151 of file OscCalculatorGeneral.cxx.

References a, ana::assert(), b, E, MECModelEnuComparisons::i, kNumFlavours, and kZeroMat().

Referenced by osc::OscCalculatorGeneral::P().

   {
     // Conversion factor for units
     E /= 4*1.267;
 
     assert(mSq.size() == kNumFlavours);
 
     ComplexMat H = kZeroMat;
 
     // Loop over rows and columns of the output
     for(unsigned int a = 0; a < kNumFlavours; ++a){
       for(unsigned int b = 0; b < kNumFlavours; ++b){
         // Form sum over mass states
         for(unsigned int i = 0; i < kNumFlavours; ++i){
           H(a, b) += U(a, i)*mSq[i]/(2*E)*std::conj(U(b, i));
         }
       }
     }
 
     return H;
   }