19 ret->Update((
unsigned char*)txt.c_str(), txt.size());
20 const int kNumParams = 8;
23 ret->Update((
unsigned char*)buf,
sizeof(
double)*kNumParams);
28 inline Eigen::Array33d
30 const Eigen::Matrix3cd & from)
32 return (to * from.adjoint()).
array().abs2();
51 __sincos(theta, &s, &c);
88 int anti = flavBefore /
abs(flavBefore);
89 int i = (
abs(flavBefore) - 12) / 2;
90 int j = (
abs(flavAfter) - 12) / 2;
91 int idx = (1-anti)/2*9 + (3 * j + i);
123 return P(flavBefore, flavAfter, E);
153 if(flavBefore < 0) anti = -1;
159 flavor_states.values,
161 flavor_states.vectors);
162 return probabilities((
abs(flavBefore)-12)/2, (
abs(flavAfter)-12)/2);
190 Eigen::Matrix3cd
const ham =
BuildHam(params);
192 for(
int anti : {-1, 1}) {
212 auto const row = Eigen::Map<const Eigen::Matrix<double, 1, 9> >(probabilities.data(), probabilities.size());
215 cache.block(
i, (1-anti)/2 * 9, 1, 9) =
row;
233 for(
int j = 0;
j < mat.cols();
j++) {
234 for(
int i = 0;
i < mat.rows();
i++) {
264 Eigen::Vector3d
const & evals,
267 const Eigen::Array3d
temp = evals.array()*
kKm2eV*(-1)*params.
L;
272 const Eigen::Array3d SIN = temp.sin();
273 const Eigen::Array3d COS = temp.cos();
274 Eigen::Vector3cd PPT;
279 std::complex<double>(COS[0], SIN[0]),
280 std::complex<double>(COS[1], SIN[1]),
281 std::complex<double>(COS[2], SIN[2]);
284 return evec * PPT.asDiagonal();
335 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3cd> eig;
336 eig.computeDirect(ham);
337 return EigenSystem{eig.eigenvectors(), eig.eigenvalues()};
347 double kr2GNe =
kK2*M_SQRT2*
kGf * params.
rho/2;
349 Eigen::Matrix3cd green_eggs(ham);
350 green_eggs.triangularView<Eigen::Upper>()/=lv;
355 green_eggs.transposeInPlace();
356 green_eggs(0,0) += kr2GNe;
359 green_eggs.adjointInPlace();
360 green_eggs(0,0) -= kr2GNe;
371 double kr2GNe =
kK2*M_SQRT2*
kGf * params.
rho/2;
373 Eigen::Matrix3cd green_eggs(ham);
374 green_eggs.triangularView<Eigen::Upper>()/=lv;
376 green_eggs.adjointInPlace();
377 green_eggs(0,0) -= kr2GNe;
389 double kr2GNe =
kK2*M_SQRT2*
kGf * params.
rho/2;
391 Eigen::Matrix3cd green_eggs(ham);
392 green_eggs.triangularView<Eigen::Upper>()/=lv;
396 green_eggs.transposeInPlace();
397 green_eggs(0,0) += kr2GNe;
412 double kr2GNe =
kK2*M_SQRT2*
kGf * params.
rho/2;
415 HLV.triangularView<Eigen::Upper>()/=lv;
418 HLV.transposeInPlace();
422 HLV.adjointInPlace();
433 double sij, cij, h00,
h11, h01;
434 std::complex<double> expCP, h02,
h12;
443 h00 = h11 * sij * sij - 1;
444 h01 = h11 * sij * cij;
445 h11 = h11 * cij * cij - 1;
449 expCP = std::complex<double>(cij, -sij);
455 h02 = (-h00 * sij * cij) * expCP;
456 h12 = (-h01 * sij) * expCP;
457 h11 -= h00 * sij * sij;
458 h00 *= cij * cij - sij * sij;
466 Eigen::Matrix3cd mHlv;
467 mHlv(0,0) = h00 - 0.5 *
h11;
468 mHlv(1,1) = 0.5 * h11 * (cij * cij - sij * sij) + 2 * real(h12) * cij * sij;
469 mHlv(2,2) = -mHlv(1,1);
471 mHlv(0,1) = h02 * sij + h01 * cij;
472 mHlv(0,2) = h02 * cij - h01 * sij;
473 mHlv(1,2) = h12 - (h11 * cij + 2 * real(h12) * sij) * sij;
TMD5 * GetParamsHash() const override
Use to check two calculators are in the same state.
Eigen::Matrix3cd BuildHam(const OscParameters ¶ms) const
Eigen::Matrix3cd AddMatterEffects(const Eigen::Matrix3cd &ham, const double &E, const int &anti, const OscParameters ¶ms) const
static const double kK2
mole/GeV^2/cm^3 to eV
Eigen::Matrix3cd AddMatterEffectsAnti(const Eigen::Matrix3cd &ham, const double &E, const OscParameters ¶ms) const
IOscCalcAdjustable * Copy() const override
void _sincos(double theta, double &s, double &c)
void PrintArrayAddresses(const double *data, int size) const
Eigen::Matrix3cd PropMatter(Eigen::Matrix3cd const &evec, Eigen::Vector3d const &evals, const OscParameters ¶ms) const
OscParameters fLastParams
void sincos(T &x, Eigen::ArrayX< U > *sx, Eigen::ArrayX< U > *cx)
double P(int flavBefore, int flavAfter, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
const XML_Char const XML_Char * data
void SetCachedEnergies(std::vector< double > const &energies)
EigenSystem Solve(Eigen::Matrix3cd const &ham) const
Eigen::Matrix3cd MatterHamiltonian(const double &E, const int &anti, const OscParameters ¶ms) const
static const double kGeV2eV
GeV to eV.
Helper struct for the cache. Might not need this.
int ChannelCacheIdx(int flavBefore, int flavAfter) const
const Eigen::Vector3d values
Eigen::Array33d Probabilities(const Eigen::Matrix3cd &to, const Eigen::Matrix3cd &from)
_OscCache< double > fCache
A re-optimized version of OscCalcPMNSOpt.
Oscillation probability calculators.
std::vector< double > energies
Array< T, Dynamic, Dynamic > probabilities
const Eigen::Matrix3cd vectors
void PrintMatrixAddresses(const double *data, int nrows, int ncols) const
_OscParameters< T > parameters
static const double kKm2eV
km to eV^-1