file src/RightHandedNeutrinos.cpp

[No description available] More…

Namespaces

Name
Gambit
TODO: see if we can use this one:
Gambit::NeutrinoBit
Eigen

Detailed Description

Author:

Date:

  • 2017 February
  • 2017 Oct
  • 2018
  • 2019
  • 2020
  • 2018
  • 2018 May

Module function definitions for NeutrinoBit exclusive for right-handed neutrinos


Authors


Source code

//   GAMBIT: Global and Modular BSM Inference Tool
//   *********************************************
///  \file
///
///  Module function definitions for NeutrinoBit 
///  exclusive for right-handed neutrinos
///
///  *********************************************
/// 
///  Authors
/// 
///  \author Suraj Krishnamurthy
///          (S.Krishnamurthy@uva.nl) 
///  \date 2017 February
///
///  \author Tomas Gonzalo
///          (t.e.gonzalo@fys.uio.no)
///  \date 2017 Oct
///  \date 2018
///  \date 2019
///  \date 2020
///
///  \author Marcin Chrzaszcz
///          (mchrzasz@cern.ch)
///  \date 2018
///
///  \author Julia Harz
///          (jharz@lpthe.jussieu.fr)
///  \date 2018 May
///
///  *********************************************

#include <cmath>
#include <iomanip>
#include <iostream>
#include <fstream>
#include <sstream>

#include "gambit/Elements/gambit_module_headers.hpp"
#include "gambit/Utils/numerical_constants.hpp"
#include "gambit/Utils/integration.hpp"
#include <unsupported/Eigen/MatrixFunctions>
#include "gambit/Utils/statistics.hpp"
#include "gambit/NeutrinoBit/NeutrinoBit_rollcall.hpp"
#include "gambit/NeutrinoBit/NeutrinoInterpolator.hpp"

//#define NEUTRINOBIT_DEBUG

using namespace Eigen;

namespace Gambit
{
  namespace NeutrinoBit
  {
    // Lambda function 
    double lambda(double a, double b, double c)
    {
      return a*a + b*b + c*c - 2*a*b - 2*b*c - 2*c*a;
    }

    // Decay with of RHNs to charged lepton and pseudoscalar meson
    // From 0705.1729, 0901.3589, 1805.08567 and 1905.00284
    double Gamma_RHN2Pplusl(double GF, double mN, double ml, double mP, double fP, double VCKM, double Usq)
    {
      double xl = pow(ml/mN,2);
      double xP = pow(mP/mN,2);

      return ( GF*GF * fP*fP * VCKM*VCKM * Usq * mN*mN*mN ) / (16*pi) * ( pow(1 - xl,2) - xP*(1 + xl) )*sqrt(lambda(1,xP,xl));
    }

    // Decay widths of RHNs, for BBN (N -> pi+ l-)
    void Gamma_RHN2piplusl(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2piplusl;
      SMInputs sminputs = *Dep::SMINPUTS;

      static const double m_pi_plus = meson_masses.pi_plus;
      static const double f_pi_plus = meson_decay_constants.pi_plus;
      // Take from the model parameters (Wolfenstein) PDG value: 0.97434
      double Vud = 1.0 - 0.5*pow(*Param["CKM_lambda"],2);
      std::vector<double> m_lep(3), gamma(3), M(3);
      // Since we scan the SLHA2 model, we take masses from it
      m_lep[0] = *Param["mE"];
      m_lep[1] = *Param["mMu"];
      m_lep[2] = *Param["mTau"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          if (M[i] > (m_pi_plus+m_lep[j]))
          {
            gamma[i] +=  Gamma_RHN2Pplusl(sminputs.GF, M[i], m_lep[j], m_pi_plus, f_pi_plus, Vud, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> K+ l-)
    void Gamma_RHN2Kplusl(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2Kplusl;
      SMInputs sminputs = *Dep::SMINPUTS;

      static const double m_K_plus = meson_masses.kaon_plus;
      static const double f_K_plus = meson_decay_constants.K_plus;
      // Take from the model parameters (Wolfenstein) PDG value: 0.22506
      double Vus = *Param["CKM_lambda"];
      std::vector<double> m_lep(3), gamma(3), M(3);
       // Since we scan the SLHA2 model, we take masses from it
      m_lep[0] = *Param["mE"];
      m_lep[1] = *Param["mMu"];
      m_lep[2] = *Param["mTau"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          if (M[i] > (m_K_plus+m_lep[j]))
          {
            gamma[i] += Gamma_RHN2Pplusl(sminputs.GF, M[i], m_lep[j], m_K_plus, f_K_plus, Vus, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> D+ l-)
    void Gamma_RHN2Dplusl(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2Dplusl;
      SMInputs sminputs = *Dep::SMINPUTS;

      static const double m_D_plus = meson_masses.D_plus;
      static const double f_D_plus = meson_decay_constants.D_plus;
      // Take from the model parameters (Wolfenstein) PDG value: 0.22492
      double Vcd = -*Param["CKM_lambda"];
      std::vector<double> m_lep(3), gamma(3), M(3);
      // Since we scan the SLHA2 model, we take masses from it
      m_lep[0] = *Param["mE"];
      m_lep[1] = *Param["mMu"];
      m_lep[2] = *Param["mTau"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          if (M[i] > (m_D_plus+m_lep[j]))
          {
            gamma[i] += Gamma_RHN2Pplusl(sminputs.GF, M[i], m_lep[j], m_D_plus, f_D_plus, Vcd, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> Ds+ l-)
    void Gamma_RHN2Dsl(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2Dsl;
      SMInputs sminputs = *Dep::SMINPUTS;

      static const double m_D_s = meson_masses.D_s;
      static const double f_D_s = meson_decay_constants.D_s;
      // Take from the model parameters (Wolfenstein) PDG value: 0.97351
      double Vcs = 1 - 0.5*pow(*Param["CKM_lambda"],2);
      std::vector<double> m_lep(3), gamma(3), M(3);
      // Since we scan the SLHA2 model, we take masses from it
      m_lep[0] = *Param["mE"];
      m_lep[1] = *Param["mMu"];
      m_lep[2] = *Param["mTau"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          if (M[i] > (m_D_s+m_lep[j]))
          {
            gamma[i] += Gamma_RHN2Pplusl(sminputs.GF, M[i], m_lep[j], m_D_s, f_D_s, Vcs, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> B+ l-)
    void Gamma_RHN2Bplusl(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2Bplusl;
      SMInputs sminputs = *Dep::SMINPUTS;

      static const double m_B_plus = meson_masses.B_plus;
      static const double f_B_plus = meson_decay_constants.B_plus;
      // Take from the model parameters (Wolfenstein) PDG value: 0.00357 (absolute value)
      double Vub = *Param["CKM_A"]*pow(*Param["CKM_lambda"],3)*sqrt(pow(*Param["CKM_rhobar"],2) + pow(*Param["CKM_etabar"],2));
      std::vector<double> m_lep(3), gamma(3), M(3);
      // Since we scan the SLHA2 model, we take masses from it
      m_lep[0] = *Param["mE"];
      m_lep[1] = *Param["mMu"];
      m_lep[2] = *Param["mTau"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          if (M[i] > (m_B_plus+m_lep[j]))
          {
            gamma[i] += Gamma_RHN2Pplusl(sminputs.GF, M[i], m_lep[j], m_B_plus, f_B_plus, Vub, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> Bc+ l-)
    void Gamma_RHN2Bcl(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2Bcl;
      SMInputs sminputs = *Dep::SMINPUTS;

      static const double m_B_c = meson_masses.B_c;
      static const double f_B_c = meson_decay_constants.B_c;
      // Take from the model parameters (Wolfenstein) PDG value: 0.0411
      double Vcb = *Param["CKM_A"]*pow(*Param["CKM_lambda"],2);
      std::vector<double> m_lep(3), gamma(3), M(3);
      // Since we scan the SLHA2 model, we take masses from it
      m_lep[0] = *Param["mE"];
      m_lep[1] = *Param["mMu"];
      m_lep[2] = *Param["mTau"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          if (M[i] > (m_B_c+m_lep[j]))
          {
            gamma[i] += Gamma_RHN2Pplusl(sminputs.GF, M[i], m_lep[j], m_B_c, f_B_c, Vcb, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay with of RHNs to neutrino and neutral pseudoscalar meson
    // From 0705.1729, 1805.08567 and 1905.00284 (mind that the latter uses Dirac neutrinos)
    double Gamma_RHN2P0nu(double GF, double mN, double mP, double fP, double Usq)
    {
      double xP = pow(mP/mN,2);

      return ( pow(GF,2) * pow(fP,2) * Usq * pow(mN,3) ) / (32*pi) * pow(1-xP,2);
    }

    // Decay widths of RHNs, for BBN (N -> pi0 nu)
    void Gamma_RHN2pi0nu(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2pi0nu;
      SMInputs sminputs = *Dep::SMINPUTS;

      static const double m_pi_0 = meson_masses.pi0;
      static const double f_pi_0 = meson_decay_constants.pi0;
      std::vector<double> gamma(3), M(3);
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        if (M[i] > m_pi_0)
        {
          for (int j=0; j<3; j++)
          {
            gamma[i] += Gamma_RHN2P0nu(sminputs.GF, M[i], m_pi_0, f_pi_0, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> eta nu)
    void Gamma_RHN2etanu(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2etanu;
      SMInputs sminputs = *Dep::SMINPUTS;

      static const double m_eta = meson_masses.eta;
      static const double f_eta = meson_decay_constants.eta;
      std::vector<double> gamma(3), M(3);
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        if (M[i] > m_eta)
        {
          for (int j=0; j<3; j++)
          {
            gamma[i] += Gamma_RHN2P0nu(sminputs.GF, M[i], m_eta, f_eta, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> eta' nu)
    void Gamma_RHN2etaprimenu(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2etaprimenu;
      SMInputs sminputs = *Dep::SMINPUTS;

      static const double m_eta_prime = meson_masses.eta_prime;
      static const double f_eta_prime = meson_decay_constants.eta_prime;
      std::vector<double> gamma(3), M(3);
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        if (M[i] > m_eta_prime)
        {
          for (int j=0; j<3; j++)
          {
            gamma[i] += Gamma_RHN2P0nu(sminputs.GF, M[i], m_eta_prime, f_eta_prime, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> etac nu)
    void Gamma_RHN2etacnu(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2etacnu;
      SMInputs sminputs = *Dep::SMINPUTS;

      static const double m_eta_c = meson_masses.eta_c;
      static const double f_eta_c = meson_decay_constants.eta_c;
      std::vector<double> gamma(3), M(3);
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        if (M[i] > m_eta_c)
        {
          for (int j=0; j<3; j++)
          {
            gamma[i] += Gamma_RHN2P0nu(sminputs.GF, M[i], m_eta_c, f_eta_c, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay with of RHNs to charged lepton and vector meson
    // From 1805.08567 and 1905.00284
    double Gamma_RHN2Vplusl(double GF, double mN, double ml, double mV, double fV, double VCKM, double Usq)
    {
      double xl = pow(ml/mN,2);
      double xV = pow(mV/mN,2);

      return ( pow(GF,2) * pow(fV,2) * pow(VCKM,2) * Usq * pow(mN,3) ) / (16*pi) * ( pow(1-xl,2) + xV*(1+xl) -2*pow(xV,2) ) * sqrt(lambda(1,xV,xl));
    }

    // Decay widths of RHNs, for BBN (N -> rho+ l-)
    void Gamma_RHN2rhoplusl(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2rhoplusl;
      SMInputs sminputs = *Dep::SMINPUTS;
 
      static const double m_rho_plus = meson_masses.rho_plus;
      static const double f_rho_plus = meson_decay_constants.rho_plus;
      // Take from the model parameters (Wolfenstein) PDG value: 0.97434
      double Vud = 1.0 - 0.5*pow(*Param["CKM_lambda"],2);
      std::vector<double> m_lep(3), gamma(3), M(3);
      // Since we scan the SLHA2 model, we take masses from it
      m_lep[0] = *Param["mE"];
      m_lep[1] = *Param["mMu"];
      m_lep[2] = *Param["mTau"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          if (M[i] > (m_rho_plus+m_lep[j]))
          {
            gamma[i] += Gamma_RHN2Vplusl(sminputs.GF, M[i], m_lep[j], m_rho_plus, f_rho_plus, Vud, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHN, for BBN (N -> D*+ l-)
    void Gamma_RHN2Dstarplusl(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2Dstarplusl;
      SMInputs sminputs = *Dep::SMINPUTS;
 
      static const double m_Dstar_plus = meson_masses.Dstar_plus;
      static const double f_Dstar_plus = meson_decay_constants.Dstar_plus;
      // Take from the model parameters (Wolfenstein) PDG value: 0.22492
      double Vcd = -*Param["CKM_lambda"];
      std::vector<double> m_lep(3), gamma(3), M(3);
      // Since we scan the SLHA2 model, we take masses from it
      m_lep[0] = *Param["mE"];
      m_lep[1] = *Param["mMu"];
      m_lep[2] = *Param["mTau"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          if (M[i] > (m_Dstar_plus+m_lep[j]))
          {
            gamma[i] += Gamma_RHN2Vplusl(sminputs.GF, M[i], m_lep[j], m_Dstar_plus, f_Dstar_plus, Vcd, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHN, for BBN (N -> Ds*+ l-)
    void Gamma_RHN2Dstarsl(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2Dstarsl;
      SMInputs sminputs = *Dep::SMINPUTS;
 
      static const double m_Dstar_s = meson_masses.Dstar_s;
      static const double f_Dstar_s = meson_decay_constants.Dstar_s;
      // Take from the model parameters (Wolfenstein) PDG value: 0.97351
      double Vcs = 1 - 0.5*pow(*Param["CKM_lambda"],2);
      std::vector<double> m_lep(3), gamma(3), M(3);
      // Since we scan the SLHA2 model, we take masses from it
      m_lep[0] = *Param["mE"];
      m_lep[1] = *Param["mMu"];
      m_lep[2] = *Param["mTau"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          if (M[i] > (m_Dstar_s+m_lep[j]))
          {
            gamma[i] += Gamma_RHN2Vplusl(sminputs.GF, M[i], m_lep[j], m_Dstar_s, f_Dstar_s, Vcs, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay with of RHNs to neutrino and neutral vector meson
    // From 1805.08567 and 1905.00284 (mind that the latter uses Dirac neutrinos)
    double Gamma_RHN2V0nu(double GF, double mN, double mV, double fV, double kappaV, double Usq)
    {
      double xV = pow(mV/mN,2);

      return ( pow(GF,2) * pow(fV,2) * pow(kappaV,2) * Usq * pow(mN,3) ) / (32*pi) * (1 + 2*xV) * pow(1 - xV,2);
    }

    // Decay widths of RHNs, for BBN (N -> rho0 nu)
    void Gamma_RHN2rho0nu(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2rho0nu;
      SMInputs sminputs = *Dep::SMINPUTS;
      double sinW2_eff = Dep::prec_sinW2_eff->central;

      static const double m_rho0 = meson_masses.rho0;
      static const double f_rho0 = meson_decay_constants.rho0;
      // kappa value from 1805.08567, it differs from that in 1905.00284
      double kappa_rho0 = 1 - 2*sinW2_eff;
      std::vector<double> gamma(3), M(3);
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        if (M[i] > m_rho0)
        {
          for (int j=0; j<3; j++)
          {
            gamma[i] += Gamma_RHN2V0nu(sminputs.GF, M[i], m_rho0, f_rho0, kappa_rho0, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> omega nu)
    void Gamma_RHN2omeganu(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2omeganu;
      SMInputs sminputs = *Dep::SMINPUTS;
      double sinW2_eff = Dep::prec_sinW2_eff->central;

      static const double m_omega = meson_masses.omega;
      static const double f_omega = meson_decay_constants.omega;
      // kappa value from 1805.08567 and 1905.00284
      double kappa_omega = (4./3.)*sinW2_eff;
      std::vector<double> gamma(3), M(3);
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        if (M[i] > m_omega)
        {
          for (int j=0; j<3; j++)
          {
            gamma[i] += Gamma_RHN2V0nu(sminputs.GF, M[i], m_omega, f_omega, kappa_omega, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> phi nu)
    void Gamma_RHN2phinu(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2phinu;
      SMInputs sminputs = *Dep::SMINPUTS;
      double sinW2_eff = Dep::prec_sinW2_eff->central;

      static const double m_phi = meson_masses.phi;
      static const double f_phi = meson_decay_constants.phi;
      // kappa value from 1805.08567 and 1905.00284
      double kappa_phi = (4./3.)*sinW2_eff - 1;
      std::vector<double> gamma(3), M(3);
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        if (M[i] > m_phi)
        {
          for (int j=0; j<3; j++)
          {
            gamma[i] += Gamma_RHN2V0nu(sminputs.GF, M[i], m_phi, f_phi, kappa_phi, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> J/psi nu)
    void Gamma_RHN2Jpsinu(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2Jpsinu;
      SMInputs sminputs = *Dep::SMINPUTS;
      double sinW2_eff = Dep::prec_sinW2_eff->central;

      static const double m_Jpsi = meson_masses.Jpsi;
      static const double f_Jpsi = meson_decay_constants.Jpsi;
      // kappa value from 1805.08567
      double kappa_Jpsi = 1 - (8./3.)*sinW2_eff;
      std::vector<double> gamma(3), M(3);
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        if (M[i] > m_Jpsi)
        {
          for (int j=0; j<3; j++)
          {
            gamma[i] += Gamma_RHN2V0nu(sminputs.GF, M[i], m_Jpsi, f_Jpsi, kappa_Jpsi, Usq(j,i));
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> nu nu nu)
    // Formula is from [arXiv:0705.1729]
    void Gamma_RHN23nu(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN23nu;
      SMInputs sminputs = *Dep::SMINPUTS;

      std::vector<double> gamma(3), M(3);
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = ( (pow(sminputs.GF,2) * pow(M[i],5)) / (192*pow(pi,3)) ) * (Usq(0,i)+Usq(1,i)+Usq(2,i));
      }
      result = gamma;
    }

    // Helper function; formula is in [arXiv:1208.4607v2]
    double S(double xa, double xb)
    {
      return sqrt((1-pow((xa+xb),2))*(1-pow((xa-xb),2)));
    }

    // Also helper function; formula is in [arXiv:1208.4607v2]
    double g(double xa, double xb)
    {
     return (1 - (7*pow(xa,2)) - (7*pow(xb,2)) - (7*pow(xa,4)) - (7*pow(xb,4)) + (12*pow(xa,2)*pow(xb,2)) - (7*pow(xa,2)*pow(xb,4)) - (7*pow(xa,4)*pow(xb,2)) + pow(xa,6) + pow(xb,6));
    }

    // Decay widths of RHNs, for BBN (N -> l+ l- nu)
    // Formula is from [arXiv:1208.4607v2]
    void Gamma_RHN2llnu(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2llnu;
      SMInputs sminputs = *Dep::SMINPUTS;

      double x_a, x_b;
      std::vector<double> m_lep(3), gamma(3), M(3);
      // Since we scan the SLHA2 model, we take masses from it
      m_lep[0] = *Param["mE"];
      m_lep[1] = *Param["mMu"];
      m_lep[2] = *Param["mTau"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          for (int k=0; k<3; k++)
          {
            if ( (j != k) and (M[i] > (m_lep[j]+m_lep[k])) )
            {
              x_a = m_lep[j]/M[i];
              x_b = m_lep[k]/M[i];
              gamma[i] += ( (pow(sminputs.GF,2)*pow(M[i],5)) / (192*pow(pi,3)) ) * Usq(j,i) * (S(x_a,x_b)*g(x_a,x_b) - (x_a < 1E-2 ? -12*pow(x_a,4) : 12*pow(x_a,4)*log( (1 - (S(x_a,x_b)*(1+pow(x_a,2)-pow(x_b,2))) - (2*pow(x_b,2)) + pow((pow(x_a,2)-pow(x_b,2)),2)) / (2*pow(x_a,2)) ) ) - (x_b < 1E-2 ? -12*pow(x_b,4) : 12*pow(x_b,4)*log( (1 - (S(x_a,x_b)*(1-pow(x_a,2)+pow(x_b,2))) - (2*pow(x_a,2)) + pow((pow(x_a,2)-pow(x_b,2)),2)) / (2*pow(x_b,2)) ) ) + (x_a < 1E-2 or x_b < 1E-2 ? -12*pow(x_a,4)*pow(x_b,4) : 12*pow(x_a,4)*pow(x_b,4)*log( (1 - (S(x_a,x_b)*(1-pow(x_a,2)-pow(x_b,2))) - (2*pow(x_a,2)) - (2*pow(x_b,2)) + pow(x_a,4) + pow(x_b,4)) / (2*pow(x_a,2)*pow(x_b,2)) ) ) );
            }
          }
        }
      }
      result = gamma;
    }

    // Helper function; formula is in [arXiv:0705.1729]
    // This function varies wrt the paper to include the x^4 factor up front and a cutoff for small x
    double L(double x)
    {
      if(x < 1E-2)
        return -pow(x,4);
      return pow(x,4)*log((1-(3*pow(x,2.0))-((1-pow(x,2.0))*sqrt(1 - (4*pow(x,2.0))))) / (pow(x,2.0)*(1+sqrt(1 - (4*pow(x,2.0))))) );
    }

    // Decay widths of RHNs, for BBN (N -> nu l+ l-)
    // Formula is from [arXiv:0705.1729]
    void Gamma_RHN2null(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2null;
      SMInputs sminputs = *Dep::SMINPUTS;
      double sinW2_eff = Dep::prec_sinW2_eff->central;

      double C1 = 0.25*(1 - (4*sinW2_eff) + (8*pow(sinW2_eff,2)));
      double C2 = 0.5*sinW2_eff*((2*sinW2_eff) - 1);
      double C3 = 0.25*(1 + (4*sinW2_eff) + (8*pow(sinW2_eff,2)));
      double C4 = 0.5*sinW2_eff*((2*sinW2_eff) + 1);
      std::vector<double> m_lep(3), gamma(3), M(3);
       // Since we scan the SLHA2 model, we take masses from it
      m_lep[0] = *Param["mE"];
      m_lep[1] = *Param["mMu"];
      m_lep[2] = *Param["mTau"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          for (int k=0; k<3; k++)
          {
            if (M[i] > (2*m_lep[k]))
            {
              double x_l = m_lep[k]/M[i];
              if (j == k)
              {
                gamma[i] += ( (pow(sminputs.GF,2)*pow(M[i],5)) / (192*pow(pi,3)) ) * Usq(j,i) * ( (C3*(((1 - (14*pow(x_l,2)) - (2*pow(x_l,4)) - (12*pow(x_l,6)))*sqrt(1 - (4*pow(x_l,2)))) + (12*(pow(x_l,4)-1)*L(x_l)))) + 4*C4*((pow(x_l,2)*(2 + (10*pow(x_l,2)) - (12*pow(x_l,4)))*sqrt(1 - (4*pow(x_l,2)))) + 6*(1.0-2*pow(x_l,2)+2*pow(x_l,4))*L(x_l)) );
              }
              else
              {
                gamma[i] += ( (pow(sminputs.GF,2)*pow(M[i],5)) / (192*pow(pi,3)) ) * Usq(j,i) * ( (C1*(((1 - (14*pow(x_l,2)) - (2*pow(x_l,4)) - (12*pow(x_l,6)))*sqrt(1 - (4*pow(x_l,2)))) + (12*(pow(x_l,4)-1)*L(x_l)))) + (4*C2*((pow(x_l,2)*(2 + (10*pow(x_l,2)) - (12*pow(x_l,4)))*sqrt(1 - (4*pow(x_l,2)))) + (6*(1.0-2*pow(x_l,2)+2*pow(x_l,4))*L(x_l)))) );
              }
            }
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> nu u ubar)
    // Formula is from 1805.08567
    void Gamma_RHN2nuuubar(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2nuuubar;
      SMInputs sminputs = *Dep::SMINPUTS;
      double sinW2_eff = Dep::prec_sinW2_eff->central;

      double C1 = 0.25*(1-(8./3)*sinW2_eff +(32./9)*pow(sinW2_eff,2));
      double C2 = (sinW2_eff/3.)*((4./3)*sinW2_eff -1);
      std::vector<double> m_uquark(3), gamma(3), M(3);
      // Since we scan the SLHA2 model, we take masses from it
      m_uquark[0] = *Param["mU"];
      m_uquark[1] = *Param["mCmC"];
      m_uquark[2] = *Param["mT"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          for (int k=0; k<3; k++)
          {
            if (M[i] > (2*m_uquark[k]))
            {
              double x = m_uquark[k]/M[i];
              // Colour factor 3
              // x^4 factor included in L function
              gamma[i] += 3*( (pow(sminputs.GF,2)*pow(M[i],5)) / (192*pow(pi,3)) ) * Usq(j,i) * (C1 * ( (1. - 14*pow(x,2) - 2*pow(x,4) -12*pow(x,6))*sqrt(1. - 4*pow(x,2)) + 12*(pow(x,4) - 1.)*L(x) ) + 4*C2 * (pow(x,2)*(2. +  10*pow(x,2) - 12*pow(x,4))*sqrt(1. - 4*pow(x,2)) + 6*(1. - 2*pow(x,2) + 2*pow(x,4))*L(x)) );
            }
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> nu d dbar)
    // Formula is from 1805.08567
    void Gamma_RHN2nuddbar(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2nuddbar;
      SMInputs sminputs = *Dep::SMINPUTS;
      double sinW2_eff = Dep::prec_sinW2_eff->central;

      double C1 = 0.25*(1-(4./3)*sinW2_eff+(8./9)*pow(sinW2_eff,2));
      double C2 = (sinW2_eff/6)*((2./3)*sinW2_eff - 1);
      std::vector<double> m_dquark(3), gamma(3), M(3);
      // Since we scan the SLHA2 model, we take masses from it
      m_dquark[0] = *Param["mD"];
      m_dquark[1] = *Param["mS"];
      m_dquark[2] = *Param["mBmB"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      // Helper function; formula is in 0705.1729 and 1805.08567
      // This function varies wrt the paper to include the x^4 factor up front and a cutoff for small x
      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          for (int k=0; k<3; k++)
          {
            if (M[i] > (2*m_dquark[k]))
            {
              double x = m_dquark[k]/M[i];
              // Colour factor 3
              // x^4 factor included in L function
              gamma[i] += 3*( (pow(sminputs.GF,2)*pow(M[i],5)) / (192*pow(pi,3)) ) * Usq(j,i) * (C1 * ( (1. - 14*pow(x,2) - 2*pow(x,4) -12*pow(x,6))*sqrt(1. - 4*pow(x,2)) + 12*(pow(x,4) - 1.)*L(x) ) + 4*C2 * (pow(x,2)*(2. +  10*pow(x,2) - 12*pow(x,4))*sqrt(1. - 4*pow(x,2)) + 6*(1. - 2*pow(x,2) + 2*pow(x,4))*L(x)) );
 
            }
          }
        }
      }
      result = gamma;
    }

    // Decay widths of RHNs, for BBN (N -> l u dbar)
    // Formula is from 1805.08567
    void Gamma_RHN2ludbar(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_RHN2ludbar;
      SMInputs sminputs = *Dep::SMINPUTS;

      // Take from the model parameters (Wolfenstein) (absolute values)
      // Values from PDG: V =  { {0.97434,0.22506,0.00357},
      //                         {0.22492,0.97351,0.0411},
      //                         {0.00875,0.0403,0.99915} };
      double Vus = *Param["CKM_lambda"];
      double Vud = 1.0 - 0.5*pow(Vus,2);
      double Vub = *Param["CKM_A"]*pow(Vus,3)*sqrt(pow(*Param["CKM_rhobar"],2) + pow(*Param["CKM_etabar"],2));
      double Vcd = Vus;
      double Vcs = Vud;
      double Vcb = *Param["CKM_A"]*pow(Vus,2);
      double Vtd = Vcb*Vus*sqrt(pow(1.0 - *Param["CKM_rhobar"],2) + pow(*Param["CKM_etabar"],2));
      double Vts = Vcb;
      double Vtb = 1;
      std::vector<std::vector<double> > V = { {Vud, Vus, Vub}, {Vcd, Vcs, Vcb}, {Vtd, Vts, Vtb} };

      std::vector<double> m_lep(3), m_uquark(3), m_dquark(3), gamma(3), M(3), decay_prod(3), two_heaviest(2);
      // Since we scan the SLHA2 model, we take masses from it
      m_lep[0] = *Param["mE"];
      m_lep[1] = *Param["mMu"];
      m_lep[2] = *Param["mTau"];
      m_uquark[0] = *Param["mU"];
      m_uquark[1] = *Param["mCmC"];
      m_uquark[2] = *Param["mT"];
      m_dquark[0] = *Param["mD"];
      m_dquark[1] = *Param["mS"];
      m_dquark[2] = *Param["mBmB"];
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2(); // |\Theta_{ij}|^2

      for (int i=0; i<3; i++)
      {
        gamma[i] = 0;
        for (int j=0; j<3; j++)
        {
          for (int k=0; k<3; k++)
          {
            for (int l=0; l<3; l++)
            {
              if (M[i] > (m_lep[j]+m_uquark[k]+m_dquark[l]))
              {
                /*decay_prod[0] = m_lep[j];
                decay_prod[1] = m_uquark[k];
                decay_prod[2] = m_dquark[l];
                two_heaviest = two_heaviest_sort(decay_prod);
                x = two_heaviest[0]/M[i];
                y = two_heaviest[1]/M[i];*/
                double xl = m_lep[j]/M[i];
                double xu = m_uquark[k]/M[i];
                double xd = m_dquark[l]/M[i];
                std::function<double(double)> I = [&xu,&xd,&xl](double x)
                {
                  return 12*(x - pow(xl,2) - pow(xd,2))*(1 + pow(xu,2) - x)/x*sqrt(lambda(x,pow(xl,2),pow(xd,2))*lambda(1,x,pow(xu,2)));
                };

                // Colour factor 3
                gamma[i] += 3*( (pow(sminputs.GF,2)*pow(M[i],5)) / (192*pow(pi,3)) ) * Usq(j,i) * pow(V[k][l],2) * Utils::integrate_cquad(I, pow(xd+xl,2), pow(1-xu,2), 0, 1e-2);
              }
            }
          }
        }
      }
      result = gamma;
    }

    // Calculates total decay width for each RHN [GeV]
    void Gamma_BBN(std::vector<double>& result)
    {
      using namespace Pipes::Gamma_BBN;

      // Charged pseudoscalar mesons
      std::vector<double> RHN2piplusl = *Dep::Gamma_RHN2piplusl;
      std::vector<double> RHN2Kplusl = *Dep::Gamma_RHN2Kplusl;
      std::vector<double> RHN2Dplusl = *Dep::Gamma_RHN2Dplusl;
      std::vector<double> RHN2Dsl = *Dep::Gamma_RHN2Dsl;
      std::vector<double> RHN2Bplusl = *Dep::Gamma_RHN2Bplusl;
      std::vector<double> RHN2Bcl = *Dep::Gamma_RHN2Bcl;

      // Neutral pseudoscalar mesons
      std::vector<double> RHN2pi0nu = *Dep::Gamma_RHN2pi0nu;
      std::vector<double> RHN2etanu = *Dep::Gamma_RHN2etanu;
      std::vector<double> RHN2etaprimenu = *Dep::Gamma_RHN2etaprimenu;
      std::vector<double> RHN2etacnu = *Dep::Gamma_RHN2etacnu;

      // Charged vector mesons
      std::vector<double> RHN2rhoplusl = *Dep::Gamma_RHN2rhoplusl;
      std::vector<double> RHN2Dstarplusl = *Dep::Gamma_RHN2Dstarplusl;
      std::vector<double> RHN2Dstarsl = *Dep::Gamma_RHN2Dstarsl;

      // Neutral vector mesons
      std::vector<double> RHN2rho0nu = *Dep::Gamma_RHN2rho0nu;
      std::vector<double> RHN2omeganu = *Dep::Gamma_RHN2omeganu;
      std::vector<double> RHN2phinu = *Dep::Gamma_RHN2phinu;
      std::vector<double> RHN2Jpsinu = *Dep::Gamma_RHN2Jpsinu;

      // Fully leptonic decays 
      std::vector<double> RHN23nu = *Dep::Gamma_RHN23nu;
      std::vector<double> RHN2llnu = *Dep::Gamma_RHN2llnu;
      std::vector<double> RHN2null = *Dep::Gamma_RHN2null;

      // Free quarks
      std::vector<double> RHN2nuuubar = *Dep::Gamma_RHN2nuuubar;
      std::vector<double> RHN2nuddbar = *Dep::Gamma_RHN2nuddbar;
      std::vector<double> RHN2ludbar = *Dep::Gamma_RHN2ludbar;

      std::vector<double> gamma_total(3), M(3);
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];
      
      // Hadronization scale 
      static const double LQCD = 1; // (GeV)

      // factor of 2 in front accounts for Majorana nature
      for (int i=0; i<3; i++)
      {
        // Fully leptonic decays are open above and below LQCD
        gamma_total[i] = 2*(RHN23nu[i] + RHN2llnu[i] + RHN2null[i]);

        // Meson decays only matter below LQCD. Beyond, only decays to free quarks matter
        if(M[i] < LQCD)
        {
          gamma_total[i] += 2*(RHN2piplusl[i] + RHN2Kplusl[i] + RHN2Dplusl[i] + RHN2Dsl[i] + RHN2Bplusl[i] + RHN2Bcl[i]);
          gamma_total[i] += 2*(RHN2pi0nu[i] + RHN2etanu[i] + RHN2etaprimenu[i] + RHN2etacnu[i]);
          gamma_total[i] += 2*(RHN2rhoplusl[i] + RHN2Dstarplusl[i] + RHN2Dstarsl[i]);
          gamma_total[i] += 2*(RHN2rho0nu[i] + RHN2omeganu[i] + RHN2phinu[i] + RHN2Jpsinu[i]);
        }
        else
        {
          gamma_total[i] += 2*(RHN2nuuubar[i] + RHN2nuddbar[i] + RHN2ludbar[i]);
        }
      }
      result = gamma_total;
    }

    // BBN constraint likelihood : lifetime must be less than 0.1s [arXiv:1202.2841]
    // This is here implemented as step function in the likelihood
    void lnL_bbn(double& result_bbn)
    {
      using namespace Pipes::lnL_bbn;
      std::vector<double> gamma = *Dep::Gamma_BBN;
      result_bbn = 0.0;
      for(int i=0; i<3; i++)
      {
        if((hbar/gamma[i])>0.1)
        {
          result_bbn = -100;
        }
      }
    }

    // Lepton universality constraint: R_{e,mu)^pi. Computation from 1502.00477. R_pi_SM from Phys. Rev. Lett 99, 231801.
    void RHN_R_pi(double& R_pi)
    {
      using namespace Pipes::RHN_R_pi;
      SMInputs sminputs = *Dep::SMINPUTS;
      static const double m_pi = meson_masses.pi_plus;
      static const double R_pi_SM = 1.2354e-4;
      double r_e_pi = pow(sminputs.mE,2)/pow(m_pi,2);
      double r_mu_pi = pow(sminputs.mMu,2)/pow(m_pi,2);
      double e_f_pi = 0.0, mu_f_pi = 0.0, d_r_pi = 1.0;
      std::vector<double> M(3), r_I_pi(3), G_e_pi = {0.0,0.0,0.0}, G_mu_pi = {0.0,0.0,0.0};
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2();

      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      for (int i=0; i<3; i++)
      {
        r_I_pi[i] = pow(M[i], 2)/pow(m_pi, 2);
 
        if(M[i] + sminputs.mMu < m_pi)
        {
          G_mu_pi[i] = ( r_mu_pi + r_I_pi[i] - pow((r_mu_pi - r_I_pi[i]), 2) ) * sqrt(1.0 - 2.0*(r_mu_pi + r_I_pi[i]) + pow(r_mu_pi - r_I_pi[i], 2)) / (r_mu_pi * pow((1.0 - r_mu_pi), 2));
        } 
        else
          G_mu_pi[i] = 0.0;

        if(M[i] + sminputs.mE < m_pi)
        {
          G_e_pi[i] = ( r_e_pi + r_I_pi[i] - pow(r_e_pi - r_I_pi[i], 2) ) * sqrt(1.0 - 2.0*(r_e_pi + r_I_pi[i]) + pow((r_e_pi - r_I_pi[i]), 2)) / (r_e_pi * pow((1.0 - r_e_pi), 2));
        }
        else
          G_e_pi[i] = 0.0;

        e_f_pi += Usq(0,i) * (G_e_pi[i] - 1.0);
        mu_f_pi += Usq(1,i) * (G_mu_pi[i] - 1.0);
  
      }

      d_r_pi = ((1.0 + e_f_pi)/(1.0 + mu_f_pi));
      R_pi = R_pi_SM * d_r_pi;
 
    }

    // Lepton universality constraint: R_(e,mu)^K. Computation from 1502.00477. R_K_SM from Phys. Rev. Lett 99, 231801.
    void RHN_R_K(double& R_K)
    {
      using namespace Pipes::RHN_R_K;
      SMInputs sminputs = *Dep::SMINPUTS;
      static const double m_K = meson_masses.kaon_plus; 
      static const double R_K_SM = 2.477e-5;
      double r_e_K = pow(sminputs.mE,2)/pow(m_K,2);
      double r_mu_K = pow(sminputs.mMu,2)/pow(m_K,2);
      double e_f_K = 0.0, mu_f_K = 0.0, d_r_K = 1.0;
      std::vector<double> M(3), r_I_K(3), G_e_K = {0.0,0.0,0.0}, G_mu_K = {0.0,0.0,0.0};
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2();

      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      for (int i=0; i<3; i++)
      {
        r_I_K[i] = pow(M[i], 2)/pow(m_K,2);

        if(M[i] + sminputs.mMu < m_K)
        {
          G_mu_K[i] = ( r_mu_K + r_I_K[i] - pow((r_mu_K - r_I_K[i]), 2) ) * sqrt(1.0 - 2.0*(r_mu_K + r_I_K[i]) + pow(r_mu_K - r_I_K[i], 2)) / (r_mu_K * pow((1.0 - r_mu_K), 2));
        } 
        else
          G_mu_K[i] = 0.0;

        if(M[i] + sminputs.mE < m_K)
        {
          G_e_K[i] = ( r_e_K + r_I_K[i] - pow((r_e_K - r_I_K[i]), 2) ) * sqrt(1.0 - 2.0*(r_e_K + r_I_K[i]) + pow((r_e_K - r_I_K[i]), 2)) / (r_e_K * pow((1.0 - r_e_K), 2));
        }
        else
          G_e_K[i] = 0.0;

        e_f_K += Usq(0,i) * (G_e_K[i] - 1.0);
        mu_f_K += Usq(1,i) * (G_mu_K[i] - 1.0);
          
      }

      d_r_K = ((1.0 + e_f_K)/(1.0 + mu_f_K));
      R_K = R_K_SM * d_r_K;
    }

    // Lepton universality constraint: R_(e,mu)^tau. Computation from 1502.00477. R_tau_SM from Int. J. Mod. Phys. A 24, 715, 2009.
    void RHN_R_tau(double& R_tau)
    {
      using namespace Pipes::RHN_R_tau;
      SMInputs sminputs = *Dep::SMINPUTS;
      static const double m_tau = sminputs.mTau;  // GeV
      static const double R_tau_SM = 0.973;
      double e_f_tau = 0.0, mu_f_tau = 0.0, d_r_tau = 1.0;
      std::vector<double> M(3);
      Matrix3d Usq = Dep::SeesawI_Theta->cwiseAbs2();

      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      for (int i=0; i<3; i++)
      {
        if(M[i] > m_tau)
        {
          e_f_tau  -= Usq(0,i);
          mu_f_tau -= Usq(1,i);
        }
        else
        {
          e_f_tau += 0.0;
          mu_f_tau += 0.0;
        }
      }

      d_r_tau = ((1.0 + mu_f_tau)/(1.0 + e_f_tau));
      R_tau = R_tau_SM * d_r_tau;
    }

    // Lepton universality from W decays
    // 0: R(W->mu nu/W->e nu) from LHCb 1608.01484
    // 1: R(W->tau nu/W->e nu) from LEP 1302.3415
    // 2: R(W->tau nu/W->mu nu) from LEP 1302.3415
    void RHN_R_W(std::vector<double> &R_W)
    {
      using namespace Pipes::RHN_R_W;
      std::vector<double> Wdecays = *Dep::W_to_l_decays;

      R_W.clear();

      R_W.push_back(Wdecays[1]/Wdecays[0]);
      R_W.push_back(Wdecays[2]/Wdecays[0]);
      R_W.push_back(Wdecays[2]/Wdecays[1]);
    }

    // Log-likelihood for the lepton universality constraint R_(e,mu)^K
    void lnL_R_K(double& result)
    {
      using namespace Pipes::lnL_R_K;
      double R_K = *Dep::R_K;
      double R_K_exp = 2.488e-5; // 1212.4012 (Phys. Lett. B 719 (2013), 326)
      double R_K_err = 0.010e-5;

      result = Stats::gaussian_loglikelihood(R_K, R_K_exp, 0.0, R_K_err, false);
    }

    // Log-likelihood for the lepton universality constraint R_(e,mu)^pi
    void lnL_R_pi(double& result)
    {
      using namespace Pipes::lnL_R_pi;
      double R_pi = *Dep::R_pi;
      double R_pi_exp = 1.2327e-4; // PDG average from Phys.Rev.Lett. 70 (1993) 17-20, Phys. Rev. Lett. 115, 071801 (2015), Phys. Rev. Lett. 68, 3000 (1992)
      double R_pi_err = 0.0023e-4;

      result = Stats::gaussian_loglikelihood(R_pi, R_pi_exp, 0.0, R_pi_err, false);
    }

    // Log-likelihood for the lepton universality constraint R_(e,mu)^tau
    void lnL_R_tau(double& result)
    {
      using namespace Pipes::lnL_R_tau;
      double R_tau = *Dep::R_tau;
      double R_tau_exp = 0.9762; // 1612.07233 (Phys. Rev. D 86, 010001) 
      double R_tau_err = 0.0028;

      result = Stats::gaussian_loglikelihood(R_tau, R_tau_exp, 0.0, R_tau_err, false);
    }
 
    // Log-likelihood for the lepton universality constraint R_(e,mu)^W 
    void lnL_R_W(double& result)
    {
      using namespace Pipes::lnL_R_W;
      std::vector<double> R_W = *Dep::R_W;

      std::vector<double> R_W_exp = {0.986, 1.043, 1.070}; // PDG 18
      std::vector<double> R_W_err = {0.013, 0.024, 0.026}; // PDG 18

      result = Stats::gaussian_loglikelihood(R_W[0], R_W_exp[0], 0.0, R_W_err[0], false);
      result += Stats::gaussian_loglikelihood(R_W[1], R_W_exp[1], 0.0, R_W_err[1], false);
      result += Stats::gaussian_loglikelihood(R_W[2], R_W_exp[2], 0.0, R_W_err[2], false);
    }

    // Calculate 0nubb half-life [1/yr] for 136Xe 0nubb detector, for the RHN model
    void RHN_Thalf_0nubb_Xe(double& result)
    {
      using namespace Pipes::RHN_Thalf_0nubb_Xe;
      double mp, A_0nubb_Xe, p2_0nubb_Xe, prefactor;
      std::vector<double> M(3);
      std::complex<double> sum = {0.0,0.0};

      // Relevant model parameters
      Matrix3cd theta = *Dep::SeesawI_Theta;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      // Nuisance parameters following the definitions in Faessler et al. 2014 (1408.6077)
      A_0nubb_Xe = runOptions->getValueOrDef<double>(8.74, "A");
      A_0nubb_Xe *= 1e-10;  // [1/yr]
      p2_0nubb_Xe = pow(runOptions->getValueOrDef<double>(183.0, "p"), 2);
      p2_0nubb_Xe *= 1e-6;  // MeV^2 --> GeV^2
      mp = 0.938;  // [GeV] (PDG 2014)

      // Lifetime equation is adopted from Faessler+14, Eq. (13)
      prefactor = A_0nubb_Xe*mp*mp/p2_0nubb_Xe/p2_0nubb_Xe;
      for (int i=0; i<3; i++)
      {
        sum+=pow(theta(0,i),2)*M[i]*p2_0nubb_Xe/(p2_0nubb_Xe+pow(M[i], 2));
      }
      result = prefactor * abs(sum) * abs(sum);
    }
    
    // Calculate 0nubb half-life [1/yr] for 76Ge 0nubb detector, for the RHN model
    void RHN_Thalf_0nubb_Ge(double& result)
    {
      using namespace Pipes::RHN_Thalf_0nubb_Ge;
      double mp, A_0nubb_Ge, p2_0nubb_Ge, prefactor;
      std::vector<double> M(3);
      std::complex<double> sum = {0.0,0.0};

      // Relevant model parameters
      Matrix3cd theta = *Dep::SeesawI_Theta;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      // Nuisance parameters following the definitions in Faessler et al. 2014 (1408.6077)
      A_0nubb_Ge = runOptions->getValueOrDef<double>(5.05, "A");
      A_0nubb_Ge *= 1e-10;  // [1/yr]
      p2_0nubb_Ge = pow(runOptions->getValueOrDef<double>(163.0, "p"), 2);
      p2_0nubb_Ge *= 1e-6;  // MeV^2 --> GeV^2
      mp = 0.938;  // [GeV] (PDG 2014)

      // Lifetime equation is adopted from Faessler+14, Eq. (13)
      prefactor = A_0nubb_Ge*mp*mp/p2_0nubb_Ge/p2_0nubb_Ge;
      for (int i=0; i<3; i++)
      {
        sum+=pow(theta(0,i),2)*M[i]*p2_0nubb_Ge/(p2_0nubb_Ge+pow(M[i], 2));
      }
      result = prefactor * abs(sum) * abs(sum);
    }

    // KamLAND-Zen: Phys. Rev. Lett 117 (2016) 082503
    void lnL_0nubb_KamLAND_Zen(double& result)
    {
      using namespace Pipes::lnL_0nubb_KamLAND_Zen;
      double tau_limit = 1.07e26;  // [yr] 90% CL

      double Thalf = *Dep::Thalf_0nubb_Xe;

      // Factor 1.28155 corresponds to one-sided UL at 90% CL
      result = Stats::gaussian_upper_limit(Thalf, 0., 0., tau_limit/1.28155, false);
    }

    // GERDA: Phys. Rev. Lett. 111 (2013) 122503
    //        Update: Nature 544 (2017) 47
    void lnL_0nubb_GERDA(double& result)
    {
      using namespace Pipes::lnL_0nubb_GERDA;
      double tau_limit = 5.3e25;  // [yr] 90% CL

      double Thalf = *Dep::Thalf_0nubb_Ge;

      // Factor 1.28155 corresponds to one-sided UL at 90% CL
      result = Stats::gaussian_upper_limit(Thalf, 0., 0., tau_limit/1.28155, false);
    }

    // Unified 0nubb likelihood
    void lnL_0nubb(double &result)
    {
      using namespace Pipes::lnL_0nubb;
      result = *Dep::lnL_0nubb_KamLAND_Zen + *Dep::lnL_0nubb_GERDA;
    }

    // Calculate mbb for 136Xe 0nubb detector, for the RHN model
    void RHN_mbb_0nubb_Xe(double& result)
    {
      using namespace Pipes::RHN_mbb_0nubb_Xe;
      double p2_0nubb_Xe;
      std::vector<double> M(3);
      std::complex<double> sum = {0.0,0.0};

      // Relevant model parameters
      Matrix3cd m_light = *Dep::m_nu;
      Matrix3cd U_light = *Dep::UPMNS;
      Matrix3cd theta = *Dep::SeesawI_Theta;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      // Nuisance parameters following the definitions in Faessler et al. 2014 (1408.6077)
      p2_0nubb_Xe = pow(runOptions->getValueOrDef<double>(178.0, "p"), 2);
      p2_0nubb_Xe *= 1e-6;  // MeV^2 --> GeV^2
      
      // mbb equation is adopted from Drewes, Eijima 2017, Eq. (14) and following
      for (int i=0; i<3; i++)
      {
        sum += pow(U_light(0,i),2)*m_light(i,i);
      }
      for (int i=0; i<3; i++)
      {
        sum += pow(theta(0,i),2)*M[i]*(p2_0nubb_Xe/(p2_0nubb_Xe + pow(M[i],2)));
      }
      result = abs(sum);
    } 
    
    // Calculate mbb for 76Ge 0nubb detector, for the RHN model
    void RHN_mbb_0nubb_Ge(double& result)
    {
      using namespace Pipes::RHN_mbb_0nubb_Ge;
      double p2_0nubb_Ge;
      std::vector<double> M(3);
      std::complex<double> sum = {0.0,0.0};

      // Relevant model parameters
      Matrix3cd m_light = *Dep::m_nu;
      Matrix3cd U_light = *Dep::UPMNS;
      Matrix3cd theta = *Dep::SeesawI_Theta;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      // Nuisance parameters following the definitions in Faessler et al. 2014 (1408.6077)
      p2_0nubb_Ge = pow(runOptions->getValueOrDef<double>(159.0, "p"), 2);
      p2_0nubb_Ge *= 1e-6;  // MeV^2 --> GeV^2
      
      // mbb equation is adopted from Drewes, Eijima 2017, Eq. (14) and following
      for (int i=0; i<3; i++)
      {
        sum += pow(U_light(0,i),2)*m_light(i,i);
      }
      for (int i=0; i<3; i++)
      {
        sum += pow(theta(0,i),2)*M[i]*(p2_0nubb_Ge/(p2_0nubb_Ge + pow(M[i],2)));
      }
      result = abs(sum);
    } 
    
    // KamLAND-Zen: Phys. Rev. Lett 117 (2016) 082503
    void lnL_mbb_0nubb_KamLAND_Zen(double& result)
    {
      using namespace Pipes::lnL_mbb_0nubb_KamLAND_Zen;
      double mbb_limit = 0.165*1e-9;  // [GeV] mbb < (0.061-0.165)eV at 90% C

      double mbb = *Dep::mbb_0nubb_Xe;

      // Factor 1.28155 corresponds to one-sided UL at 90% CL
      result = Stats::gaussian_upper_limit(mbb, 0., 0., mbb_limit/1.28155, false);
    }

    // GERDA: Phys. Rev. Lett. 111 (2013) 122503
    //        Update: Nature 544 (2017) 47
    void lnL_mbb_0nubb_GERDA(double& result)
    {
      using namespace Pipes::lnL_mbb_0nubb_GERDA;
      double mbb_limit = 0.33*1e-9;  // [GeV] mbb < (0.15-0.33)eV at 90% CL

      double mbb = *Dep::mbb_0nubb_Ge;

      // Factor 1.28155 corresponds to one-sided UL at 90% CL
      result = Stats::gaussian_upper_limit(mbb, 0., 0., mbb_limit/1.28155, false);
    }

    // Unified 0nubb likelihood based on mbb
    void lnL_mbb_0nubb(double &result)
    {
      using namespace Pipes::lnL_mbb_0nubb;
      result = *Dep::lnL_mbb_0nubb_KamLAND_Zen + *Dep::lnL_mbb_0nubb_GERDA;
    }

    // Helper function to fill all needed experimental info for CKM unitarity
    void fill_ckm_exp(Matrix3cd &Theta, double GF, triplet<double> (&Vus_exp)[8], triplet<double> &Vud_exp, double (&f)[8])
    {
      // Vus
      // Experimental values determined for K and tau decays. From table 1 in 1502.00477
      double V_us_exp[] = {0.2163, 0.2166, 0.2155, 0.2160, 0.2158, 0.2262, 0.2214, 0.2173};
      double err_V_us_exp[] = {0.0006, 0.0006, 0.0013, 0.0011, 0.0014, 0.0013, 0.0022, 0.0022};
      double f_plus = 0.959;
      double err_f_plus = 0.005;
      for(int i=0; i<5; i++)
      {
        V_us_exp[i] /= f_plus;
        err_V_us_exp[i] = sqrt(pow(err_V_us_exp[i] / f_plus,2) + pow(V_us_exp[i] * err_f_plus / f_plus, 2));
      }
 
      // Fill Vus_exp triplet
      for(int i=0; i<8; i++)
      {
        Vus_exp[i].central = V_us_exp[i];
        Vus_exp[i].upper = err_V_us_exp[i];
        Vus_exp[i].lower = err_V_us_exp[i];
      }

      // Vud
      // Combined value from the PDG
      static const double V_ud_exp = 0.97417;
      static const double err_V_ud_exp = 0.00021;
 
      // Fill Vud_exp triplet
      Vud_exp.central = V_ud_exp;
      Vud_exp.upper = err_V_ud_exp;
      Vud_exp.lower = err_V_ud_exp;

      // f
      Matrix3d ThetaNorm = (Theta * Theta.adjoint()).real();
      double Gmu = sqrt(GF*GF*(1. - ThetaNorm(1,1) - ThetaNorm(0,0)));

      f[0] = pow(GF/Gmu,2)*(1 - ThetaNorm(0,0));
      f[1] = f[0];
      f[2] = f[0];
      f[3] = pow(GF/Gmu,2)*(1 - ThetaNorm(1,1));
      f[4] = f[3];
      f[5] = 1 + ThetaNorm(1,1);
      f[6] = 1 + ThetaNorm(0,0) + ThetaNorm(1,1) - ThetaNorm(2,2);
      f[7] = 1 + 0.2*ThetaNorm(0,0) - 0.9*ThetaNorm(1,1) - 0.2*ThetaNorm(2,2);
    }
     
    // CKM unitarity constraint: optimize on Vus from Theta [PDG 2016]
    void calc_Vus(double& result_Vus)
    {
      using namespace Pipes::calc_Vus;
      SMInputs sminputs = *Dep::SMINPUTS;
      Matrix3cd Theta = *Dep::SeesawI_Theta;

      // Fill experimental data
      triplet<double> V_us_exp[8];
      triplet<double> V_ud_exp;
      double f[8];
      fill_ckm_exp(Theta, sminputs.GF, V_us_exp, V_ud_exp, f);
      
      // For the minimalization it's much better to transform the Vud experimental result to Vus the same as we do for theory and minimize only Vus
      double V_us_from_Vud=sqrt(1.-V_ud_exp.central*V_ud_exp.central);
      double err_V_us_from_Vud_exp= ( (V_ud_exp.central)/(sqrt(1-V_ud_exp.central*V_ud_exp.central))  ) * V_ud_exp.upper;
                                             
      double est_Vus = 0.;
      double sum_numerator = 0;
      double sum_denominator = 0;
      for (int i=0; i<7; i++)
      {
        sum_numerator += (V_us_exp[i].central/f[i]) / (V_us_exp[i].upper*V_us_exp[i].upper /( f[i]*f[i]));
        sum_denominator += 1./(V_us_exp[i].upper*V_us_exp[i].upper /( f[i]*f[i]));
      }
      // now vud
      sum_numerator += (V_us_from_Vud /f[0]) / ( (err_V_us_from_Vud_exp*err_V_us_from_Vud_exp)/( f[0]*f[0]));
      sum_denominator +=  1./( (err_V_us_from_Vud_exp*err_V_us_from_Vud_exp) /( f[0]*f[0]) );

      est_Vus = sum_numerator/sum_denominator;
      
      result_Vus = est_Vus;
    }

    // CKM unitarity constraint: V_ud should lie within 3sigma of the world average [PDG 2016]
    void lnL_ckm_Vusmin(double& result_ckm)
    {
      using namespace Pipes::lnL_ckm_Vusmin;
      SMInputs sminputs = *Dep::SMINPUTS;
      Matrix3cd Theta = *Dep::SeesawI_Theta;
      double V_us = *Dep::calc_Vus;

      // Fill experimental data
      triplet<double> V_us_exp[8];
      triplet<double> V_ud_exp;
      double f[8];
      fill_ckm_exp(Theta, sminputs.GF, V_us_exp, V_ud_exp, f);
 
      double chi2 = 0.0;
      for (int i=0; i<7; i++)
        chi2 += pow( (sqrt(pow(V_us,2)*f[i]) - V_us_exp[i].central) / V_us_exp[i].upper, 2);
      // According to 1407.6607 the correction for Vud is the same as K->pi e nu (f[0])
      double V_ub = 3.94e-3;
      chi2 += pow( (sqrt((1 - pow(V_us,2) - pow(V_ub,2))*f[0]) - V_ud_exp.central)/ V_ud_exp.upper/1.20, 2);
      result_ckm = -0.5*chi2;
    }
  
    // CKM unitarity constraint: V_ud should lie within 3sigma of the world average [PDG 2016]
    void lnL_ckm_Vus(double& result_ckm)
    {
      using namespace Pipes::lnL_ckm_Vus;
      SMInputs sminputs = *Dep::SMINPUTS;
      Matrix3cd Theta = *Dep::SeesawI_Theta;
      double V_us = *Param["CKM_lambda"];

      // Fill experimental data
      triplet<double> V_us_exp[8];
      triplet<double> V_ud_exp;
      double f[8];
      fill_ckm_exp(Theta, sminputs.GF, V_us_exp, V_ud_exp, f);

      double chi2 = 0.0;
      for (int i=0; i<7; i++)
        chi2 += pow( (sqrt(pow(V_us,2)*f[i]) - V_us_exp[i].central) / V_us_exp[i].upper, 2);
      // According to 1407.6607 the correction for Vud is the same as K->pi e nu (f[0])
      chi2 += pow( (sqrt((1 - pow(V_us,2))*f[0]) - V_ud_exp.central)/ V_ud_exp.upper, 2);
      result_ckm = -0.5*chi2;
    }
    
    // Likelihood contribution from PIENU; searched for extra peaks in the spectrum of pi -> mu + nu. Constrains |U_ei|^2 at 90% in the mass range 60-129 MeV. [Phys. Rev. D, 84 (052002), 2011] (arXiv:1106.4055)
    void lnL_pienu(double& result)
    {
      using namespace Pipes::lnL_pienu;

      // Mass range of experiment
      static const double low_lim = 0.06060759493670887;  // GeV
      static const double upp_lim = 0.12926582278481014;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(3);

      mixing_sq[0] = *Dep::Ue1;
      mixing_sq[1] = *Dep::Ue2;
      mixing_sq[2] = *Dep::Ue3;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/pienu.csv");

      // Assume Gaussian errors with zero mean and that limits scale as |U|^2.
      result = 0;
      for(int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
        {
          result += -0.5*log(2.0*pi/1.28/1.28);
        }
        else
        {
          U[i] = s.eval(M[i]);
          result += Stats::gaussian_upper_limit(mixing_sq[i]/U[i], 0, 0, 1/1.28, false);  // exp_error = abs(exp_value - 90CL_value), exp_value = 0, 1.28: 90% CL limit for half-Gaussian.
        }
      }
    }

    // Likelihood contribution from PS191, electron sector; looked for charged tracks originating from RHN decays: nu_r -> l(-) + l(+) + nu / l + pi / e + pi(+) + pi(0). Constrains |U_ei|^2 at 90% in the mass range 20-450 MeV. Function also incorporates a later re-interpretation of the data to account for neutral current interaction (ignored in original) as well as the RHNs' Majorana nature. [Original: Phys. Lett. B, 203(3):332-334, 1988][Re-interp.: JHEP, 2012(6):1-27 (arXiv:1112.3319)]
    void lnL_ps191_e(double& result)
    {
      using namespace Pipes::lnL_ps191_e;

      // Mass range of experiment
      static const double low_lim = 0.011810586;  // GeV
      static const double upp_lim = 0.4491907842;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(3);
      double c_e = 0.5711;
      double c_mu = 0.1265;
      double c_tau = 0.1265;

      mixing_sq[0] = *Dep::Ue1 * ((c_e * *Dep::Ue1) + (c_mu * *Dep::Um1) + (c_tau * *Dep::Ut1));
      mixing_sq[1] = *Dep::Ue2 * ((c_e * *Dep::Ue2) + (c_mu * *Dep::Um2) + (c_tau * *Dep::Ut2));
      mixing_sq[2] = *Dep::Ue3 * ((c_e * *Dep::Ue3) + (c_mu * *Dep::Um3) + (c_tau * *Dep::Ut3));
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/ps191_e.csv");

      // Assume scaling with |U|^4, zero bkg, number of events at 90% CL is
      // reverse engineered.  We assume that lnL = mu_sig is a faithful
      // approximation to the true Poisson likelihood.
      result = 0;
      for (int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
          result += 0;
        else
        {
          U[i] = s.eval(M[i]);
          result += -2.44*(mixing_sq[i]/pow(U[i],2));
        }
      }
    }

    // Likelihood contribution from PS191, muon sector. Constrains |U_(mu,i)|^2 at 90% in the mass range 20-450 MeV. Description & references above.
    void lnL_ps191_mu(double& result)
    {
      using namespace Pipes::lnL_ps191_mu;

      // Mass range of experiment
      static const double low_lim = 0.0103348894;  // GeV
      static const double upp_lim = 0.3610401587;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(3);
      double c_e = 0.5711;
      double c_mu = 0.1265;
      double c_tau = 0.1265;

      mixing_sq[0] = *Dep::Um1 * ((c_e * *Dep::Ue1) + (c_mu * *Dep::Um1) + (c_tau * *Dep::Ut1));
      mixing_sq[1] = *Dep::Um2 * ((c_e * *Dep::Ue2) + (c_mu * *Dep::Um2) + (c_tau * *Dep::Ut2));
      mixing_sq[2] = *Dep::Um3 * ((c_e * *Dep::Ue3) + (c_mu * *Dep::Um3) + (c_tau * *Dep::Ut3));
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/ps191_mu.csv");

      // Assume scaling with |U|^4, zero bkg, number of events at 90% CL is
      // reverse engineered.  We assume that lnL = mu_sig is a faithful
      // approximation to the true Poisson likelihood.
      result = 0;
      for (int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
          result += 0;
        else
        {
          U[i] = s.eval(M[i]);
          result += -2.44*(mixing_sq[i]/pow(U[i],2));
        }
      }
    }

    // Likelihood contribution from CHARM, electron sector; searched for charged and neutral current decays of RHNs. Constrains |U_ei|^2 at 90% in the mass range 0.5-2.8 GeV. [Phys. Lett. B, 166(4):473-478, 1986]
    void lnL_charm_e(double& result)
    {
      using namespace Pipes::lnL_charm_e;

      // Mass range of experiment
      static const double low_lim = 0.1595725833606568;  // GeV
      static const double upp_lim = 2.0814785578102417;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(3);
      double c_e = 0.5711;
      double c_mu = 0.1265;
      double c_tau = 0.1265;

      mixing_sq[0] = *Dep::Ue1 * ((c_e * *Dep::Ue1) + (c_mu * *Dep::Um1) + (c_tau * *Dep::Ut1));
      mixing_sq[1] = *Dep::Ue2 * ((c_e * *Dep::Ue2) + (c_mu * *Dep::Um2) + (c_tau * *Dep::Ut2));
      mixing_sq[2] = *Dep::Ue3 * ((c_e * *Dep::Ue3) + (c_mu * *Dep::Um3) + (c_tau * *Dep::Ut3));
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/charm_e.csv");

      // Assume scaling with |U|^4, zero bkg, number of events at 90% CL is
      // reverse engineered.  We assume that lnL = mu_sig is a faithful
      // approximation to the true Poisson likelihood.
      result = 0;
      for (int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
          result += 0;
        else
        {
          U[i] = s.eval(M[i])/sqrt(2);  // Division by sqrt(2) to account for Majorana nature.
          result += -2.44*(mixing_sq[i]/pow(U[i],2));
        }
      }
    }

    // Likelihood contribution from CHARM, muon sector. Constrains |U_(mu,i)|^2 at 90% in the mass range 0.5-2.8 GeV. Description & references above.
    void lnL_charm_mu(double& result)
    {
      using namespace Pipes::lnL_charm_mu;

      // Mass range of experiment
      static const double low_lim = 0.4483153374997989;  // GeV
      static const double upp_lim = 1.9231171483448785;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(3);
      double c_e = 0.5711;
      double c_mu = 0.1265;
      double c_tau = 0.1265;

      mixing_sq[0] = *Dep::Um1 * ((c_e * *Dep::Ue1) + (c_mu * *Dep::Um1) + (c_tau * *Dep::Ut1));
      mixing_sq[1] = *Dep::Um2 * ((c_e * *Dep::Ue2) + (c_mu * *Dep::Um2) + (c_tau * *Dep::Ut2));
      mixing_sq[2] = *Dep::Um3 * ((c_e * *Dep::Ue3) + (c_mu * *Dep::Um3) + (c_tau * *Dep::Ut3));
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/charm_mu.csv");

      // Assume scaling with |U|^4, zero bkg, number of events at 90% CL is
      // reverse engineered.  We assume that lnL = mu_sig is a faithful
      // approximation to the true Poisson likelihood.
      result = 0;
      for (int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
          result += 0;
        else
        {
          U[i] = s.eval(M[i])/sqrt(2);  // Division by sqrt(2) to account for Majorana nature.
          result += -2.44*(mixing_sq[i]/pow(U[i],2));
        }
      }
    }

    // Likelihood contribution from DELPHI's short-lived RHN analysis; searched for charged and neutral current decays of RHNs. Constrains |U_ei|^2, |U_(mu,i)|^2 as well as |U_(tau,i)|^2 at 95% in the mass range 3.5-50 GeV. [Z. Phys. C, 74(1):57-71, 1997]
    void lnL_delphi_short_lived(double& result)
    {
      using namespace Pipes::lnL_delphi_short_lived;

      // Mass range of experiment
      static const double low_lim = 1.8102188251700203;  // GeV
      static const double upp_lim = 80.00000000000006;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(9);

      mixing_sq[0] = *Dep::Ue1;  // This is |U_{e1}|^2 etc
      mixing_sq[1] = *Dep::Ue2;
      mixing_sq[2] = *Dep::Ue3;
      mixing_sq[3] = *Dep::Um1;
      mixing_sq[4] = *Dep::Um2;
      mixing_sq[5] = *Dep::Um3;
      mixing_sq[6] = *Dep::Ut1;
      mixing_sq[7] = *Dep::Ut2;
      mixing_sq[8] = *Dep::Ut3;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/delphi_short_lived.csv");

      // Assume scaling with |U|^4, zero bkg, number of events at 95% CL is
      // reverse engineered.  We assume that lnL = mu_sig is a faithful
      // approximation to the true Poisson likelihood.
      result = 0;
      for (int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
          result += 0;
        else
        {
          for (int j=i; j<9; j+=3)
          {
            U[i] = s.eval(M[i]);
            result += -3.09*(pow(mixing_sq[j]/U[i],2));
          }
        }
      }
    }

    // Likelihood contribution from DELPHI's long-lived RHN analysis; searched for charged and neutral current decays of RHNs. Constrains |U_ei|^2, |U_(mu,i)|^2 as well as |U_(tau,i)|^2 at 95% in the mass range 0.5-4.2 GeV. [Z. Phys. C, 74(1):57-71, 1997]
    void lnL_delphi_long_lived(double& result)
    {
      using namespace Pipes::lnL_delphi_long_lived;

      // Mass range of experiment
      static const double low_lim = 0.4383563;  // GeV
      static const double upp_lim = 4.1954595;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(9);

      mixing_sq[0] = *Dep::Ue1;  // This is |U_{e1}|^2 etc
      mixing_sq[1] = *Dep::Ue2;
      mixing_sq[2] = *Dep::Ue3;
      mixing_sq[3] = *Dep::Um1;
      mixing_sq[4] = *Dep::Um2;
      mixing_sq[5] = *Dep::Um3;
      mixing_sq[6] = *Dep::Ut1;
      mixing_sq[7] = *Dep::Ut2;
      mixing_sq[8] = *Dep::Ut3;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/delphi_long_lived.csv");

      // Assume scaling with |U|^4, zero bkg, number of events at 95% CL is
      // reverse engineered.  We assume that lnL = mu_sig is a faithful
      // approximation to the true Poisson likelihood.
      result = 0;
      for (int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
          result += 0;
        else
        {
          for (int j=i; j<9; j+=3)
          {
            U[i] = s.eval(M[i]);
            result += -3.09*(pow(mixing_sq[j]/U[i],2));
          }
        }
      }
    }

    // Likelihood contribution from ATLAS, electron sector; looked at the production and decay chain: pp -> W*(+-) -> l(+-) + nu_r. nu_r then decays into an on-shell W and a lepton; the W decays primarily into a qq pair. Constrains |U_ei|^2 at 95% in the mass range 50-500 GeV. [JHEP, 07:162, 2015 (arXiv:1506.06020)]
    void lnL_atlas_e(double& result)
    {
      using namespace Pipes::lnL_atlas_e;

      // Mass range of experiment
      static const double low_lim = 100.1041668;  // GeV
      static const double upp_lim = 476.1458333;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(3);

      mixing_sq[0] = *Dep::Ue1;
      mixing_sq[1] = *Dep::Ue2;
      mixing_sq[2] = *Dep::Ue3;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/atlas_e.csv");

      // Assume Gaussian errors with zero mean and that limits scale as |U|^4.
      result = 0;
      for(int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
        {
          result += -0.5*log(2.0*pi/1.64/1.64);
        }
        else
        {
          U[i] = s.eval(M[i]);
          result += Stats::gaussian_upper_limit(pow(mixing_sq[i]/U[i],2), 0, 0, 1/1.64, false);  // exp_error = abs(exp_value - 95CL_value), exp_value = 0, 1.64: 95% CL limit for half-Gaussian.
        }
      }
    }

    // Likelihood contribution from ATLAS, muon sector. Constrains |U_(mu,i)|^2 at 95% in the mass range 50-500 GeV. Description & references above.
    void lnL_atlas_mu(double& result)
    {
      using namespace Pipes::lnL_atlas_mu;

      // Mass range of experiment
      static const double low_lim = 101.89094090419824;  // GeV
      static const double upp_lim = 500.76903192219294;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(3);

      mixing_sq[0] = *Dep::Um1;
      mixing_sq[1] = *Dep::Um2;
      mixing_sq[2] = *Dep::Um3;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/atlas_mu.csv");

      // Assume Gaussian errors with zero mean and that limits scale as |U|^4.
      result = 0;
      for(int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
        {
          result += -0.5*log(2.0*pi/1.64/1.64);
        }
        else
        {
          U[i] = s.eval(M[i]);
          result += Stats::gaussian_upper_limit(pow(mixing_sq[i]/U[i],2), 0, 0, 1/1.64, false); // exp_error = abs(exp_value - 95CL_value), exp_value = 0, 1.64: 95% CL limit for half-Gaussian.
        }
      }
    }

    // Likelihood contribution from E949; used the kaon decay: K(+) -> mu(+) + nu_r. Constrains |U_(mu,i)|^2 at 90% in the mass range 175-300 MeV. [Phys. Rev. D, 91, 052001 (2015) (arXiv:1411.3963v2)]
    void lnL_e949(double& result)
    {
      using namespace Pipes::lnL_e949;

      // Mass range of experiment
      static const double low_lim = 0.1794613032227713;  // GeV
      static const double upp_lim = 0.2995365796283227;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(3);

      mixing_sq[0] = *Dep::Um1;
      mixing_sq[1] = *Dep::Um2;
      mixing_sq[2] = *Dep::Um3;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/e949.csv");

      // Assume Gaussian errors with zero mean and that limits scale as |U|^2.
      result = 0;
      for(int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
        {
          result += -0.5*log(2.0*pi/1.28/1.28);
        }
        else
        {
          U[i] = s.eval(M[i])/sqrt(2);  // Division by sqrt(2) to account for Majorana nature.
          result += Stats::gaussian_upper_limit(mixing_sq[i]/U[i], 0, 0,  1/1.28, false); // exp_error = abs(exp_value - 90CL_value), exp_value = 0, 1.28: 90% CL limit for half-Gaussian.
        }
      }
    }

    // Likelihood contribution from NuTeV; used RHN decays into muonic final states (mu + mu + nu / mu + e + nu / mu + pi / mu + rho). Constrains |U_(mu,i)|^2 at 90% CL in the mass range 0.25-2 GeV. [Phys. Rev. Lett., 83:4943-4946, 1999 (arXiv:hep-ex/9908011)]
    void lnL_nutev(double& result)
    {
      using namespace Pipes::lnL_nutev;

      // Mass range of experiment
      static const double low_lim = 0.2116390354;  // GeV
      static const double upp_lim = 2.0161957132;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(3);

      mixing_sq[0] = *Dep::Um1;
      mixing_sq[1] = *Dep::Um2;
      mixing_sq[2] = *Dep::Um3;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/nutev.csv");

      // Assume scaling with |U|^4, zero bkg, number of events at 90% CL is
      // reverse engineered.  We assume that lnL = mu_sig is a faithful
      // approximation to the true Poisson likelihood.
      result = 0;
      for (int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
          result += 0;
        else
        {
          U[i] = s.eval(M[i]);
          result += -2.44*(pow(mixing_sq[i]/U[i],2));
        }
      }
    }

    // Likelihood contribution from a re-interpretation of CHARM data; assumes tau mixing is dominant. Constrains |U_(tau,i)|^2 at 90% CL in the mass range 10-290 MeV. [Phys. Lett. B, 550(1-2):8-15, 2002 (arXiv:hep-ph/0208075)]
    void lnL_charm_tau(double& result)
    {
      using namespace Pipes::lnL_charm_tau;

      // Mass range of experiment
      static const double low_lim = 0.010685579196217497;  // GeV
      static const double upp_lim = 0.288745725563687;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(3);

      mixing_sq[0] = *Dep::Ut1;
      mixing_sq[1] = *Dep::Ut2;
      mixing_sq[2] = *Dep::Ut3;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/tau.csv");

      // Assume Gaussian errors with zero mean and that limits scale as |U|^4.
      result = 0;
      for(int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
        {
          result += -0.5*log(2*pi/1.28/1.28);
        }
        else
        {
          U[i] = s.eval(M[i])/sqrt(2);  // Division by sqrt(2) to account for Majorana nature.
          result += Stats::gaussian_upper_limit(pow(mixing_sq[i]/U[i],2), 0, 0, 1/1.28, false); // exp_error = abs(exp_value - 95CL_value), exp_value = 0, 1.28: 90% CL limit for half-Gaussian.
        }
      }
    }

    // Likelihood contribution from LHC (CMS), electron sector; looked at the production and decay chain: pp -> W*(+-) -> l(+-) + nu_r. nu_r then decays into an on-shell W and a lepton; the W decays primarily into a qq pair. Constrains |U_ei|^2 at 95% in the mass range 1-1.2e3 GeV. [arXiv:1802.02965v1]
    void lnL_lhc_e(double& result)
    {
      using namespace Pipes::lnL_lhc_e;

      // Mass range of experiment
      static const double low_lim = 1.0293246;  // GeV
      static const double upp_lim = 1e3;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(3);

      mixing_sq[0] = *Dep::Ue1;
      mixing_sq[1] = *Dep::Ue2;
      mixing_sq[2] = *Dep::Ue3;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/lhc_e.csv");

      // Assume Gaussian errors with zero mean and that limits scale as |U|^4.
      result = 0;
      for(int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
        {
          result += -0.5*log(2.0*pi/1.64/1.64);
        }
        else
        {
          U[i] = s.eval(M[i]);
          result += Stats::gaussian_upper_limit(pow(mixing_sq[i]/U[i],2), 0, 0, 1/1.64, false);  // exp_error = abs(exp_value - 95CL_value), exp_value = 0, 1.64: 95% CL limit for half-Gaussian.
        }
      }
    }

    // Likelihood contribution from LHC (CMS), muon sector. Constrains |U_(mu,i)|^2 at 95% in the mass range 1-1.2e3 GeV. Description and references above.
    void lnL_lhc_mu(double& result)
    {
      using namespace Pipes::lnL_lhc_mu;

      // Mass range of experiment
      static const double low_lim = 1.0145564;  // GeV
      static const double upp_lim = 985.652549;  // GeV
      std::vector<double> M(3), U(3), mixing_sq(3);

      mixing_sq[0] = *Dep::Um1;
      mixing_sq[1] = *Dep::Um2;
      mixing_sq[2] = *Dep::Um3;
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      static NeutrinoInterpolator s("NeutrinoBit/data/lhc_mu.csv");

      // Assume Gaussian errors with zero mean and that limits scale as |U|^4.
      result = 0;
      for(int i=0; i<3; i++)
      {
        if ( (M[i] < low_lim) or (M[i] > upp_lim) )
        {
          result += -0.5*log(2.0*pi/1.64/1.64);
        }
        else
        {
          U[i] = s.eval(M[i]);
          result += Stats::gaussian_upper_limit(pow(mixing_sq[i]/U[i],2), 0, 0, 1/1.64, false);  // exp_error = abs(exp_value - 95CL_value), exp_value = 0, 1.64: 95% CL limit for half-Gaussian.
        }
      }
    }

    // Squared matrix element |Theta(1,1)|^2
    void Ue1(double& Ue1_sq)
    {
      using namespace Pipes::Ue1;
      Ue1_sq = ((*Dep::SeesawI_Theta).cwiseAbs2())(0,0);

      #ifdef NEUTRINOBIT_DEBUG
        cout << "Ue1 = " << Ue1_sq << endl;
      #endif

      double upper_limit = runOptions->getValueOrDef<double>(-1, "upper_limit");
      double lower_limit = runOptions->getValueOrDef<double>(-1, "lower_limit");

      if( (upper_limit != -1 and Ue1_sq > upper_limit) or (lower_limit != -1 and Ue1_sq < lower_limit) )
      {
        std::ostringstream msg;
        msg << "Coupling outside of given limits";
        logger() << msg.str() << EOM;
        invalid_point().raise(msg.str());
      }
    }

    // Phase of matrix element Theta(1,1)
    void Ue1_phase(double& Ue1_p)
    {
      using namespace Pipes::Ue1_phase;
      Ue1_p = std::arg((*Dep::SeesawI_Theta)(0,0));
    }

    // Squared matrix element |Theta(2,1)|^2
    void Um1(double& Um1_sq)
    {
      using namespace Pipes::Um1;
      Um1_sq = (Dep::SeesawI_Theta->cwiseAbs2())(1,0);

      #ifdef NEUTRINOBIT_DEBUG
        cout << "Um1 = " << Um1_sq << endl;
      #endif

      double upper_limit = runOptions->getValueOrDef<double>(-1, "upper_limit");
      double lower_limit = runOptions->getValueOrDef<double>(-1, "lower_limit");

      if( (upper_limit != -1 and Um1_sq > upper_limit) or (lower_limit != -1 and Um1_sq < lower_limit) )
      {
        std::ostringstream msg;
        msg << "Coupling outside of given limits";
        logger() << msg.str() << EOM;
        invalid_point().raise(msg.str());
      }
 
    }

    // Phase of matrix element Theta(2,1)
    void Um1_phase(double& Um1_p)
    {
      using namespace Pipes::Um1_phase;
      Um1_p = std::arg((*Dep::SeesawI_Theta)(1,0));
    }

    // Squared matrix element |Theta(3,1)|^2
    void Ut1(double& Ut1_sq)
    {
      using namespace Pipes::Ut1;
      Ut1_sq = (Dep::SeesawI_Theta->cwiseAbs2())(2,0);

      #ifdef NEUTRINOBIT_DEBUG
        cout << "Ut1 = " << Ut1_sq << endl;
      #endif

      double upper_limit = runOptions->getValueOrDef<double>(-1, "upper_limit");
      double lower_limit = runOptions->getValueOrDef<double>(-1, "lower_limit");

      if( (upper_limit != -1 and Ut1_sq > upper_limit) or (lower_limit != -1 and Ut1_sq < lower_limit) )
      {
        std::ostringstream msg;
        msg << "Coupling outside of given limits";
        logger() << msg.str() << EOM;
        invalid_point().raise(msg.str());
      }
 
    }

    // Phase of matrix element Theta(3,1)
    void Ut1_phase(double& Ut1_p)
    {
      using namespace Pipes::Ut1_phase;
      Ut1_p = std::arg((*Dep::SeesawI_Theta)(2,0));
    }

    // Squared matrix element |Theta(1,2)|^2
    void Ue2(double& Ue2_sq)
    {
      using namespace Pipes::Ue2;
      Ue2_sq = (Dep::SeesawI_Theta->cwiseAbs2())(0,1);

      #ifdef NEUTRINOBIT_DEBUG
        cout << "Ue2 = " << Ue2_sq << endl;
      #endif

      double upper_limit = runOptions->getValueOrDef<double>(-1, "upper_limit");
      double lower_limit = runOptions->getValueOrDef<double>(-1, "lower_limit");

      if( (upper_limit != -1 and Ue2_sq > upper_limit) or (lower_limit != -1 and Ue2_sq < lower_limit) )
      {
        std::ostringstream msg;
        msg << "Coupling outside of given limits";
        logger() << msg.str() << EOM;
        invalid_point().raise(msg.str());
      }
 
    }

    // Phase of matrix element Theta(1,2)
    void Ue2_phase(double& Ue2_p)
    {
      using namespace Pipes::Ue2_phase;
      Ue2_p = std::arg((*Dep::SeesawI_Theta)(0,1));
    }

    // Squared matrix element |Theta(2,2)|^2
    void Um2(double& Um2_sq)
    {
      using namespace Pipes::Um2;
      Um2_sq = (Dep::SeesawI_Theta->cwiseAbs2())(1,1);
 
      #ifdef NEUTRINOBIT_DEBUG
        cout << "Um2 = " << Um2_sq << endl;
      #endif

      double upper_limit = runOptions->getValueOrDef<double>(-1, "upper_limit");
      double lower_limit = runOptions->getValueOrDef<double>(-1, "lower_limit");

      if( (upper_limit != -1 and Um2_sq > upper_limit) or (lower_limit != -1 and Um2_sq < lower_limit) )
      {
        std::ostringstream msg;
        msg << "Coupling outside of given limits";
        logger() << msg.str() << EOM;
        invalid_point().raise(msg.str());
      }
 
    }

    // Phase of matrix element Theta(2,2)
    void Um2_phase(double& Um2_p)
    {
      using namespace Pipes::Um2_phase;
      Um2_p = std::arg((*Dep::SeesawI_Theta)(1,1));
    }

    // Squared matrix element |Theta(3,2)|^2
    void Ut2(double& Ut2_sq)
    {
      using namespace Pipes::Ut2;
      Ut2_sq = (Dep::SeesawI_Theta->cwiseAbs2())(2,1);

      #ifdef NEUTRINOBIT_DEBUG
        cout << "Ut2 = " << Ut2_sq << endl;
      #endif

      double upper_limit = runOptions->getValueOrDef<double>(-1, "upper_limit");
      double lower_limit = runOptions->getValueOrDef<double>(-1, "lower_limit");

      if( (upper_limit != -1 and Ut2_sq > upper_limit) or (lower_limit != -1 and Ut2_sq < lower_limit) )
      {
        std::ostringstream msg;
        msg << "Coupling outside of given limits";
        logger() << msg.str() << EOM;
        invalid_point().raise(msg.str());
      }
 
    }

    // Phase of matrix element Theta(3,2)
    void Ut2_phase(double& Ut2_p)
    {
      using namespace Pipes::Ut2_phase;
      Ut2_p = std::arg((*Dep::SeesawI_Theta)(2,1));
    }

    // Squared matrix element |Theta(1,3)|^2
    void Ue3(double& Ue3_sq)
    {
      using namespace Pipes::Ue3;
      Ue3_sq = (Dep::SeesawI_Theta->cwiseAbs2())(0,2);

      #ifdef NEUTRINOBIT_DEBUG
        cout << "Ue3 = " << Ue3_sq << endl;
      #endif

      double upper_limit = runOptions->getValueOrDef<double>(-1, "upper_limit");
      double lower_limit = runOptions->getValueOrDef<double>(-1, "lower_limit");

      if( (upper_limit != -1 and Ue3_sq > upper_limit) or (lower_limit != -1 and Ue3_sq < lower_limit) )
      {
        std::ostringstream msg;
        msg << "Coupling outside of given limits";
        logger() << msg.str() << EOM;
        invalid_point().raise(msg.str());
      }
 
    }

    // Phase of matrix element Theta(1,3)
    void Ue3_phase(double& Ue3_p)
    {
      using namespace Pipes::Ue3_phase;
      Ue3_p = std::arg((*Dep::SeesawI_Theta)(0,2));
    }

    // Squared matrix element |Theta(2,3)|^2
    void Um3(double& Um3_sq)
    {
      using namespace Pipes::Um3;
      Um3_sq = (Dep::SeesawI_Theta->cwiseAbs2())(1,2);

      #ifdef NEUTRINOBIT_DEBUG
        cout << "Um3 = " << Um3_sq << endl;
      #endif

      double upper_limit = runOptions->getValueOrDef<double>(-1, "upper_limit");
      double lower_limit = runOptions->getValueOrDef<double>(-1, "lower_limit");

      if( (upper_limit != -1 and Um3_sq > upper_limit) or (lower_limit != -1 and Um3_sq < lower_limit) )
      {
        std::ostringstream msg;
        msg << "Coupling outside of given limits";
        logger() << msg.str() << EOM;
        invalid_point().raise(msg.str());
      }
 
    }

    // Phase of matrix element Theta(2,3)
    void Um3_phase(double& Um3_p)
    {
      using namespace Pipes::Um3_phase;
      Um3_p = std::arg((*Dep::SeesawI_Theta)(1,2));
    }

    // Squared matrix element |Theta(3,3)|^2
    void Ut3(double& Ut3_sq)
    {
      using namespace Pipes::Ut3;
      Ut3_sq = (Dep::SeesawI_Theta->cwiseAbs2())(2,2);

      #ifdef NEUTRINOBIT_DEBUG
        cout << "Ut3 = " << Ut3_sq << endl;
      #endif

      double upper_limit = runOptions->getValueOrDef<double>(-1, "upper_limit");
      double lower_limit = runOptions->getValueOrDef<double>(-1, "lower_limit");

      if( (upper_limit != -1 and Ut3_sq > upper_limit) or (lower_limit != -1 and Ut3_sq < lower_limit) )
      {
        std::ostringstream msg;
        msg << "Coupling outside of given limits";
        logger() << msg.str() << EOM;
        invalid_point().raise(msg.str());
      }
 
    }

    // Phase of matrix element Theta(3,3)
    void Ut3_phase(double& Ut3_p)
    {
      using namespace Pipes::Ut3_phase;
      Ut3_p = std::arg((*Dep::SeesawI_Theta)(2,2));
    }

    // Step-function log likelihood on the perturbativity of the Yukawas
    void perturbativity_likelihood(double &lnL)
    {
      using namespace Pipes::perturbativity_likelihood;
      SMInputs sminputs = *Dep::SMINPUTS;

      Matrix3d MN;
      MN << *Param["M_1"], 0, 0,
            0, *Param["M_2"], 0,
            0, 0, *Param["M_3"];

      double vev= 1. / sqrt(sqrt(2.)*sminputs.GF);

      // Yukawa coupling |F|^2 from eq 26 in 1502.00477
      Matrix3cd F2 = 1.0/pow(vev,2) * *Dep::SeesawI_Theta * Dep::SeesawI_Theta->adjoint() * MN * MN;
      
      lnL = 0;
      for(int i=0; i<3; i++)
        for(int j=0; j<3; j++)
          if(F2(i,j).real() >= 4*pi) 
          {
            std::ostringstream msg;
            msg << "Yukawas not perturbative; point invalidated.";
            logger() << msg.str() << EOM;
            invalid_point().raise(msg.str());
            return ;
          }
    }

    // Artificial slide likelihood on couplings and masses
    void coupling_slide(double &lnL)
    {
      using namespace Pipes::coupling_slide;
      int I = runOptions->getValueOrDef<int>(1, "I");
      int flavour = runOptions->getValueOrDef<int>(1, "i");
      double slope = runOptions->getValueOrDef<double>(1, "slope");
      double mslope = runOptions->getValueOrDef<double>(0, "mslope");
      std::vector<double> M(3);
      M[0] = *Param["M_1"];
      M[1] = *Param["M_2"];
      M[2] = *Param["M_3"];

      if (flavour > 0)
      {
        double U = (Dep::SeesawI_Theta->cwiseAbs2())(flavour-1,I-1);
        lnL = slope*log10(std::min(U, 1.)) + mslope*log10(M[I-1]);
      }
      else
      {
        double U1 = (Dep::SeesawI_Theta->cwiseAbs2())(0, I-1);
        double U2 = (Dep::SeesawI_Theta->cwiseAbs2())(1, I-1);
        double U3 = (Dep::SeesawI_Theta->cwiseAbs2())(2, I-1);
        lnL = 0;
        lnL += slope*log10(std::min(U1, 1.)) + mslope*log10(M[I-1]);
        lnL += slope*log10(std::min(U2, 1.)) + mslope*log10(M[I-1]);
        lnL += slope*log10(std::min(U3, 1.)) + mslope*log10(M[I-1]);
      }
    }
  }
}

Updated on 2024-07-18 at 13:53:33 +0000