Add three flavour barger propagator (#139)

* to be used in testing for comparison with tensor based propagator

Author: ewanwm

python/binding.cpp

@@ -419,4 +419,39 @@ void initTesting(py::module &m)
             py::arg("energy"), py::arg("i"), py::arg("j")
         )
     ;
+
+    py::class_<testing::ThreeFlavourBarger>(m_testing, "ThreeFlavourBarger")
+        .def(py::init<>())
+        .def("set_params", &testing::ThreeFlavourBarger::setParams, 
+            py::arg("m1"), py::arg("m2"), py::arg("m3"), py::arg("theta12"), py::arg("theta13"), py::arg("theta23"), py::arg("deltaCP"), py::arg("baseline"), py::arg("density") = (float)-999.9, py::arg("anti_neutrino") = false
+        )
+        .def("alpha", &testing::ThreeFlavourBarger::alpha,
+            "Calculates alpha term used in calculating the mass eigenvalues",
+            py::arg("energy")
+        )
+        .def("beta", &testing::ThreeFlavourBarger::beta,
+            "Calculates beta term used in calculating the mass eigenvalues",
+            py::arg("energy")
+        )
+        .def("gamma", &testing::ThreeFlavourBarger::gamma,
+            "Calculates gamma term used in calculating the mass eigenvalues",
+            py::arg("energy")
+        )
+        .def("calculate_effective_m2", &testing::ThreeFlavourBarger::calculateEffectiveM2,
+            "Calculates the effective hamiltonian eigenvalues (the M^2) in matter",
+            py::arg("energy"), py::arg("index")
+        )
+        .def("get_hamiltonian_element", &testing::ThreeFlavourBarger::getHamiltonianElement,
+            "Calculates an element of the Hamiltonian",
+            py::arg("energy"), py::arg("a"), py::arg("b")
+        )
+        .def("get_transition_matrix_element", &testing::ThreeFlavourBarger::getTransitionMatrixElement,
+            "Calculates an element of the transition matrix from one mass eigenstate to another due to the presense of matter",
+            py::arg("energy"), py::arg("a"), py::arg("b")
+        )
+        .def("calculate_prob", &testing::ThreeFlavourBarger::calculateProb,
+            "Calculate probability of transitioning from state i to state j for a given energy",
+            py::arg("energy"), py::arg("i"), py::arg("j")
+        )
+    ;
 }

tests/barger-propagator.hpp

@@ -1,6 +1,9 @@
 #pragma once
 
+#include <array>
+
 #include <cmath>
+#include <complex>
 
 #include <iostream>
 #include <nuTens/propagator/constants.hpp>
@@ -195,6 +198,249 @@ class TwoFlavourBarger
     bool _antiNeutrino;
 };
 
+class ThreeFlavourBarger 
+{
+  public:
+
+    // set the parameters of this propagator
+    // negative density values will be interpreted as propagating in vacuum
+    inline void setParams(
+        double m1, double m2, double m3, 
+        double theta12, double theta13, double theta23, 
+        double deltaCP,
+        double baseline, double density = -999.9f, bool antiNeutrino = false)
+    {
+        _m1 = m1;
+        _m2 = m2;
+        _m3 = m3;
+        _theta12 = theta12;
+        _theta13 = theta13;
+        _theta23 = theta23;
+        _deltaCP = deltaCP;
+        _baseline = baseline;
+        _density = density;
+        _antiNeutrino = antiNeutrino;
+
+        // fill the mass array
+        masses[0] = _m1;
+        masses[1] = _m2;
+        masses[2] = _m3;
+
+        // fill the PMNS matrix elements
+        pmnsMatrix[0][0] = std::complex<double>(std::cos(theta12) * std::cos(theta13), 0.0);
+        pmnsMatrix[0][1] = std::complex<double>(std::sin(theta12) * std::cos(theta13), 0.0);
+        pmnsMatrix[0][2] = std::sin(theta13) * std::exp(std::complex<double>(0.0, -1.0) * deltaCP);
+
+        pmnsMatrix[1][0] = -std::sin(theta12) * std::cos(theta23) - std::cos(theta12) * std::sin(theta23) * std::sin(theta13) * std::exp(std::complex<double>(0.0, 1.0) * deltaCP);
+        pmnsMatrix[1][1] = std::cos(theta12) * std::cos(theta23) - std::sin(theta12) * std::sin(theta23) * std::sin(theta13) * std::exp(std::complex<double>(0.0, 1.0) * deltaCP);
+        pmnsMatrix[1][2] = std::complex<double>(std::sin(theta23) * std::cos(theta13), 0.0);
+
+        pmnsMatrix[2][0] = std::sin(theta12) * std::sin(theta23) - std::cos(theta12) * std::cos(theta23) * std::sin(theta13) * std::exp(std::complex<double>(0.0, 1.0) * deltaCP);
+        pmnsMatrix[2][1] = -std::cos(theta12) * std::sin(theta23) - std::sin(theta12) * std::cos(theta23) * std::sin(theta13) * std::exp(std::complex<double>(0.0, 1.0) * deltaCP);
+        pmnsMatrix[2][2] = std::complex<double>(std::cos(theta23) * std::cos(theta13), 0.0);
+
+    };
+
+    /// calculate the alpha factor used in the eigenvalue computation
+    [[nodiscard]] inline double alpha(float energy) const 
+    {
+        float dmsq12 = _m1 * _m1 - _m2 * _m2;
+        float dmsq13 = _m1 * _m1 - _m3 * _m3;
+
+        double ret;
+        
+        if (_antiNeutrino) {
+            ret = - 2.0 * constants::Groot2 * energy * _density + dmsq12 + dmsq13;
+        }
+        else {
+            ret = 2.0 * constants::Groot2 * energy * _density + dmsq12 + dmsq13;
+        }
+
+        return ret;
+    }
+
+    /// calculate the beta factor used in the eigenvalue computation
+    [[nodiscard]] inline double beta(float energy) const 
+    {
+        double dmsq12 = _m1 * _m1 - _m2 * _m2;
+        double dmsq13 = _m1 * _m1 - _m3 * _m3;
+
+        double ret;
+
+        if (_antiNeutrino) {
+            ret = (
+                dmsq12 * dmsq13 + 
+                - 2.0 * constants::Groot2 * energy * _density * (
+                    dmsq12 * (1.0 - std::abs(pmnsMatrix[0][1]) * std::abs(pmnsMatrix[0][1])) +
+                    dmsq13 * (1.0 - std::abs(pmnsMatrix[0][2]) * std::abs(pmnsMatrix[0][2]))
+                )
+            );
+        }
+        else {
+            ret = (
+                dmsq12 * dmsq13 + 
+                2.0 * constants::Groot2 * energy * _density * (
+                    dmsq12 * (1.0 - std::abs(pmnsMatrix[0][1]) * std::abs(pmnsMatrix[0][1])) +
+                    dmsq13 * (1.0 - std::abs(pmnsMatrix[0][2]) * std::abs(pmnsMatrix[0][2]))
+                )
+            );
+        }
+
+        return ret;
+    }
+
+    /// calculate the gamma factor used in the eigenvalue computation
+    [[nodiscard]] inline double gamma(float energy) const 
+    {
+        float dmsq12 = _m1 * _m1 - _m2 * _m2;
+        float dmsq13 = _m1 * _m1 - _m3 * _m3;
+
+        double ret;
+        if (_antiNeutrino) {
+            ret = -2 * constants::Groot2 * energy * _density * dmsq12 * dmsq13 * std::abs(pmnsMatrix[0][0]) * std::abs(pmnsMatrix[0][0]);
+        }
+        else {
+            ret = 2 * constants::Groot2 * energy * _density * dmsq12 * dmsq13 * std::abs(pmnsMatrix[0][0]) * std::abs(pmnsMatrix[0][0]);
+        }
+
+        return ret;
+    }
+
+    /// calculate effective M^2 values (eigenvalues of the hamiltonian) due to matter effects
+    /// @param energy The neutrino energy
+    /// @param index The index of the eigenvalue. should be [0-2]
+    [[nodiscard]] inline double calculateEffectiveM2(float energy, int index) const 
+    {
+        float a = alpha(energy);
+        float b = beta(energy);
+        float c = gamma(energy);
+
+        // calculate argument of arccos
+        float arg = (2.0 * a*a*a - 9.0 * a*b + 27.0 * c) / ( 2.0 * std::pow( a*a - 3.0 * b, 3.0 / 2.0) ); 
+
+        // calculate the coefficient of the cos term
+        float coeff = - (2.0 / 3.0) * std::sqrt( a*a - 3.0 * b );
+
+        return coeff * std::cos( ( 1.0 / 3.0 ) * ( std::acos(arg) + index * 2.0 * M_PI ) ) + _m1 * _m1 - a / 3.0;
+    }
+
+    /// @brief Calculate an element of the hamiltonian
+    /// @param energy The neutrino energy
+    /// @param k Row
+    /// @param j Column
+    /// @return Matrix element
+    [[nodiscard]] inline std::complex<double> getHamiltonianElement(float energy, int a, int b) const {
+        
+        std::complex<double> ret = 0.0;
+
+        if ( a == b ) {
+            ret += masses[a] * masses[a] / (2.0 * energy);
+        }
+
+        if (_antiNeutrino) {
+            ret += constants::Groot2 * _density * std::conj(pmnsMatrix[0][b]) * pmnsMatrix[0][a];
+        }
+        else {
+            ret -= constants::Groot2 * _density * pmnsMatrix[0][b] * std::conj(pmnsMatrix[0][a]);
+        }
+
+        return ret;
+    }
+
+    /// @brief Calculate an element of the "X" transition matrix matrix (equation 11 in Barger et al)
+    /// @param energy The neutrino energy
+    /// @param a Row
+    /// @param b Column
+    /// @return Matrix element
+    [[nodiscard]] inline std::complex<double> getTransitionMatrixElement(float energy, int a, int b) const {
+
+        std::complex<double> ret = 0.0;
+
+        for (int k = 0; k < 3; k++) {
+
+            std::complex<double> numerator = 4.0 * energy * energy * (
+                getHamiltonianElement(energy, a, 0) * getHamiltonianElement(energy, 0, b) +
+                getHamiltonianElement(energy, a, 1) * getHamiltonianElement(energy, 1, b) +
+                getHamiltonianElement(energy, a, 2) * getHamiltonianElement(energy, 2, b) 
+            );
+
+            std::complex<double> constant = 1.0;
+            std::complex<double> denominator = 1.0;
+
+            for (int j = 0; j < 3; j++) {
+
+                if (j == k) continue;
+
+                numerator -= 2.0 * energy * getHamiltonianElement(energy, a, b) * calculateEffectiveM2(energy, j);
+                denominator *= calculateEffectiveM2(energy, k) - calculateEffectiveM2(energy, j);
+                constant *= calculateEffectiveM2(energy, j);
+
+            }
+
+            std::complex<double> prod = numerator;
+
+            if (a == b) 
+                prod += constant;
+        
+            prod /= denominator;
+
+            std::complex<double> exponential = std::exp(-std::complex<double>(0.0, 1.0) * calculateEffectiveM2(energy, k) * _baseline * 2.0 * M_PI / (2.0 * energy));
+
+            ret += prod * exponential;
+        }
+
+        return ret;
+    }
+
+    /// @brief calculate oscillation probability from flavour alpha to flavour beta
+    /// @param energy neutrino energy
+    /// @param alpha initial flavour index
+    /// @param beta final flavour index
+    /// @return oscillation probability
+    [[nodiscard]] inline double calculateProb(float energy, int alpha, int beta) const {
+
+        std::complex<double> ret = 0.0;
+
+        for (int i = 0; i < 3; i++) {
+
+            for (int j = 0; j < 3; j++) {
+            
+                if (_antiNeutrino) {
+                    ret += std::conj(pmnsMatrix[beta][i] * getTransitionMatrixElement(energy, i, j)) * pmnsMatrix[alpha][j];
+                }
+                else {
+                    ret += pmnsMatrix[alpha][i] * getTransitionMatrixElement(energy, i, j) * std::conj(pmnsMatrix[beta][j]);
+                }
+
+            }
+        
+        }
+
+        return std::abs(ret) * std::abs(ret);
+    }
+    
+
+  private:
+    // oscillation parameters
+    double _m1;
+    double _m2;
+    double _m3;
+    double _theta12;
+    double _theta13;
+    double _theta23;
+    double _deltaCP;
+
+    // other parameters
+    double _baseline;
+    double _density;
+
+    // anti-neutrino flag
+    bool _antiNeutrino;
+
+    std::array<std::array<std::complex<double>, 3>, 3> pmnsMatrix;
+    std::array<double, 3> masses;
+};
+
 } // namespace testing
 
 } // namespace nuTens