Merge branch 'master' into gqd

MainPage.md

@@ -15,6 +15,7 @@ Available classes are:
 - **[PMNS_Deco](@ref OscProb::PMNS_Deco):** Oscillations with 3 flavours including a simple decoherence model
 - **[PMNS_LIV](@ref OscProb::PMNS_LIV):** Oscillations with 3 flavours including Lorentz Invariance Violations
 - **[PMNS_Decay](@ref OscProb::PMNS_Decay):** Oscillations with 3 flavours including neutrino decays
+- **[PMNS_NUNM](@ref OscProb::PMNS_NUNM):** Oscillations with 3 flavours including non-unitary neutrino mixing
 - **[Absorption](@ref OscProb::Absorption):** Computes absorption probabilities for high-energy neutrinos
 
 A few example macros on how to use OscProb are available in a [tutorial](https://github.com/joaoabcoelho/OscProb/tree/master/tutorial) directory.

README.md

@@ -17,6 +17,7 @@ Available classes are:
 - **PMNS_Deco:** Oscillations with 3 flavours including a simple decoherence model
 - **PMNS_LIV:** Oscillations with 3 flavours including Lorentz Invariance Violations
 - **PMNS_Decay:** Oscillations with 3 flavours including neutrino decays
+- **PMNS_NUNM:** Oscillations with 3 flavours including non-unitary neutrino mixing
 - **Absorption:** Computes absorption probabilities for high-energy neutrinos
 
 A few example macros on how to use OscProb are available in a [tutorial](tutorial) directory.

inc/PMNS_NUNM.h

@@ -0,0 +1,83 @@
+///////////////////////////////////////////////////////////////////////////////
+/// \class OscProb::PMNS_NUNM
+///
+/// \brief Implementation of oscillations of neutrinos in matter in a
+///        three-neutrino framework with Non unitary Neutrino Mixing (NUNM).
+///
+/// This class expands the PMNS_Fast class including a general NU mixing matrix
+///
+///
+/// The non unitarity effect is parametrized by dimensionless quantities
+/// alpha which quantify the deviation from unitarity with respect to the
+/// standard mixing
+///
+/// Reference: https://arxiv.org/pdf/2111.00329.pdf
+///
+/// \sa PMNS_Fast
+///
+/// \author [email protected]
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef PMNS_NUNM_H
+#define PMNS_NUNM_H
+
+#include "Definitions.h"
+#include "PMNS_Fast.h"
+#include <Eigen/Dense>
+
+namespace OscProb {
+
+  class PMNS_NUNM : public PMNS_Fast {
+    public:
+      PMNS_NUNM(int scale = 0); ///< Constructor
+      virtual ~PMNS_NUNM();     ///< Destructor
+
+      /// Set any given NUNM parameter
+      virtual void SetAlpha(int i, int j, double val,
+                            double phase); // i, j between 0 and 2
+
+      /// Get any given NUNM parameter
+      virtual complexD GetAlpha(int i, int j);
+
+      /// Set the NUNM parameters all at once
+      void SetNUNM(double alpha_11, double alpha_21, double alpha_31,
+                   double alpha_22, double alpha_32, double alpha_33);
+
+      // Set the diagonal real NUNM parameters
+      virtual void SetAlpha_11(double a); ///< Set alpha_11 parameter
+      virtual void SetAlpha_22(double a); ///< Set alpha_22 parameter
+      virtual void SetAlpha_33(double a); ///< Set alpha_33 parameter
+
+      // Set the off-diagonal complex NUNM parameters
+      virtual void    SetAlpha_21(double a,
+                                  double phi); ///< Set alpha_21 parameter
+      virtual void    SetAlpha_31(double a,
+                                  double phi); ///< Set alpha_31 parameter
+      virtual void    SetAlpha_32(double a,
+                                  double phi); ///< Set alpha_32 parameter
+      virtual void    SetFracVnc(double f);
+      virtual matrixD ProbMatrix(int nflvi, int nflvf);
+
+    protected:
+      virtual void UpdateHam();
+      virtual void Propagate();
+      virtual void PropagatePath(NuPath p);
+      vectorC      ApplyAlphaDagger(vectorC fState);
+      vectorC      ApplyAlpha(vectorC fState);
+      void         InitMatrix();
+
+      int     fscale;
+      double  fracVnc; // set fraction of matter potential affecting NC
+      vectorC fNuStateBuffer;
+
+      Eigen::Matrix<std::complex<double>, 3, 3> X;
+      Eigen::Matrix<std::complex<double>, 3, 3> Alpha;
+      Eigen::Matrix<std::complex<double>, 3, 3> V;
+      Eigen::Matrix<std::complex<double>, 3, 3> Ham;
+  };
+
+} // namespace OscProb
+
+#endif
+
+////////////////////////////////////////////////////////////////////////