Commit fcedd5d9 authored by Grégoire Uhlrich's avatar Grégoire Uhlrich
Browse files

A_tau fixed in pMSSM

parent 1291c2c3
......@@ -250,24 +250,24 @@ void print_libdiagonalization_cppdata(std::ostream &out) {
out << " std::vector<complex> m_dagger_m = dot(m_data_herm, m_data);\n";
out << "\n";
out << " std::vector<real> mass2 = m_mass;\n";
out << " diagonalize(m_m_dagger, m_transfer, m_mass, true);\n";
out << " diagonalize(m_dagger_m, m_transfer2, mass2, true);\n";
out << " sort(m_transfer, m_mass);\n";
out << " sort(m_transfer2, mass2);\n";
out << " diagonalize(m_m_dagger, m_transfer2, m_mass, true);\n";
out << " diagonalize(m_dagger_m, m_transfer, mass2, true);\n";
out << " sort(m_transfer2, m_mass);\n";
out << " sort(m_transfer, mass2);\n";
out << "\n";
out << " for (size_t i = 0; i != m_N; ++i) {\n";
out << " assert(ABS(m_mass[i] - mass2[i]) / ((m_mass[i] == 0) ? 1 : m_mass[i])\n";
out << " < precision);\n";
out << " < precision);\n";
out << " m_mass[i] = SQRT(m_mass[i]);\n";
out << " }\n";
out << "\n";
out << " auto diagMass = dot(dot(hermitian(m_transfer), m_data), m_transfer2);\n";
out << " auto diagMass = dot(dot(hermitian(m_transfer2), m_data), m_transfer);\n";
out << " for (size_t i = 0; i != m_N; ++i) \n";
out << " if (REAL(diagMass[index(i, i)]*m_mass[i]) < 0) \n";
out << " for (size_t j = 0; j != m_N; ++j)\n";
out << " m_transfer2[index(j, i)] *= -1;\n";
out << " m_transfer[index(j, i)] *= -1;\n";
out << "\n";
out << " for (auto &el : m_transfer)\n";
out << " for (auto &el : m_transfer2)\n";
out << " el = CONJ(el);\n";
out << "}\n";
out << "\n";
......
// This file is part of MARTY.
//
// MARTY is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// MARTY is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with MARTY. If not, see <https://www.gnu.org/licenses/>.
#include <cassert>
#include "libdiagonalization.h"
#include <gsl/gsl_complex.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex_math.h>
#include <iostream>
#ifdef QUAD
#include <quadmath.h>
#define ABS fabsq
#define SQRT sqrtq
#define REAL(arg) crealq(arg)
#define IMAG(arg) cimagq(arg)
#define CONJ(arg) conjq(arg)
#else
#define ABS std::abs
#define SQRT std::sqrt
#define REAL(arg) (arg).real()
#define IMAG(arg) (arg).imag()
#define CONJ(arg) std::conj(arg)
#endif
Diagonalizer::real Diagonalizer::precision = 1e-5;
void Diagonalizer::applyDiagonalization(
std::initializer_list<complex> massMatrix,
std::initializer_list<complex*> transfer,
std::initializer_list<real*> masses
)
{
Diagonalizer diag(massMatrix);
std::vector<complex> const &t = diag.transfer();
size_t i = 0;
for (auto iter = transfer.begin(); iter != transfer.end(); ++iter)
**iter = t[i++];
i = 0;
std::vector<real> const &m = diag.mass();
for (auto iter = masses.begin(); iter != masses.end(); ++iter)
**iter = m[i++];
}
void Diagonalizer::applyBiDiagonalization(
std::initializer_list<complex> massMatrix,
std::initializer_list<complex*> transfer,
std::initializer_list<complex*> transfer2,
std::initializer_list<real*> masses
)
{
Diagonalizer diag(massMatrix, true);
std::vector<complex> const &t = diag.transfer();
size_t i = 0;
for (auto iter = transfer.begin(); iter != transfer.end(); ++iter)
**iter = t[i++];
i = 0;
std::vector<complex> const &t2 = diag.transfer2();
for (auto iter = transfer2.begin(); iter != transfer2.end(); ++iter)
**iter = t2[i++];
i = 0;
std::vector<real> const &m = diag.mass();
for (auto iter = masses.begin(); iter != masses.end(); ++iter)
**iter = m[i++];
}
Diagonalizer::Diagonalizer(
std::initializer_list<complex> massMatrix,
bool t_bidiag
)
:m_data(massMatrix),
m_N(std::sqrt(m_data.size())),
m_bidiag(t_bidiag),
m_computed(false)
{
assert(m_N * m_N == massMatrix.size());
}
std::vector<Diagonalizer::complex> const &Diagonalizer::transfer()
{
updateDiagonalization();
return m_transfer;
}
std::vector<Diagonalizer::complex> const &Diagonalizer::transfer2()
{
updateDiagonalization();
assert(m_bidiag);
return m_transfer2;
}
std::vector<Diagonalizer::real> const &
Diagonalizer::mass()
{
updateDiagonalization();
return m_mass;
}
void Diagonalizer::updateDiagonalization()
{
if (!m_computed) {
(m_bidiag) ? bidiagonalize() : diagonalize();
m_computed = true;
}
}
void Diagonalizer::diagonalize(
std::vector<complex> massMatrix,
std::vector<complex> &transfer,
std::vector<real> &mass,
bool bidiagonalization
)
{
transfer = std::vector<complex>(m_N*m_N, 0);
mass = std::vector<real>(m_N, 0);
gsl_eigen_hermv_workspace* workspace = gsl_eigen_hermv_alloc(m_N);
if (!bidiagonalization) {
// Symmetrize the mass matrix (if it has not been done yet)
const size_t N = std::sqrt(massMatrix.size());
for (size_t i = 0; i != N; ++i)
for (size_t j = i+1; j < N; ++j) {
massMatrix[index(i, j)] =
(massMatrix[index(i, j)] + massMatrix[index(j, i)]) / 2.;
massMatrix[index(j, i)] = massMatrix[index(i, j)];
}
}
gsl_matrix_complex* init_gsl = getGSLMassMatrix(massMatrix);
gsl_vector* eigenValues = gsl_vector_alloc(m_N);
gsl_matrix_complex* transfer_gsl = gsl_matrix_complex_alloc(m_N, m_N);
// Diagonalization
gsl_eigen_hermv(init_gsl, eigenValues, transfer_gsl, workspace);
// Gathering results
loadGSLMatrix(transfer_gsl, transfer);
loadGSLVector(eigenValues, mass);
// Freeing allocated memory by gsl
gsl_matrix_complex_free(transfer_gsl);
gsl_vector_free (eigenValues);
gsl_matrix_complex_free(init_gsl);
gsl_eigen_hermv_free (workspace);
}
void Diagonalizer::swap(
const size_t i,
const size_t j,
std::vector<complex> &transfer,
std::vector<real> &mass
)
{
const size_t N = mass.size();
for (size_t k = 0; k != N; ++k) {
std::swap(transfer[index(k, i)], transfer[index(k, j)]);
}
}
void Diagonalizer::sort(
std::vector<complex> &transfer,
std::vector<real> &mass
)
{
for (size_t i = 0; i != mass.size(); ++i) {
size_t i_min = i;
for (size_t j = i+1; j < mass.size(); ++j) {
if (ABS(mass[j]) < ABS(mass[i_min]))
i_min = j;
}
if (i != i_min) {
swap(i, i_min, transfer, mass);
std::swap(mass[i_min], mass[i]);
}
}
}
void Diagonalizer::sort(
std::vector<complex> &transfer,
std::vector<complex> &transfer2,
std::vector<real> &mass
)
{
for (size_t i = 0; i != mass.size(); ++i) {
size_t i_min = i;
for (size_t j = i+1; j < mass.size(); ++j) {
if (ABS(mass[j]) < ABS(mass[i_min]))
i_min = j;
}
if (i != i_min) {
swap(i, i_min, transfer, mass);
swap(i, i_min, transfer2, mass);
std::swap(mass[i_min], mass[i]);
}
}
}
void Diagonalizer::diagonalize()
{
diagonalize(m_data, m_transfer, m_mass);
sort(m_transfer, m_mass);
}
void Diagonalizer::bidiagonalize()
{
std::vector<complex> m_data_herm = hermitian(m_data);
std::vector<complex> m_m_dagger = dot(m_data, m_data_herm);
std::vector<complex> m_dagger_m = dot(m_data_herm, m_data);
std::vector<real> mass2 = m_mass;
diagonalize(m_m_dagger, m_transfer, m_mass, true);
diagonalize(m_dagger_m, m_transfer2, mass2, true);
sort(m_transfer, m_mass);
sort(m_transfer2, mass2);
for (size_t i = 0; i != m_N; ++i) {
assert(ABS(m_mass[i] - mass2[i]) / ((m_mass[i] == 0) ? 1 : m_mass[i])
< precision);
m_mass[i] = SQRT(m_mass[i]);
}
for (auto &el : m_transfer)
el = CONJ(el);
sort(m_transfer, m_transfer2, m_mass);
}
gsl_matrix_complex *Diagonalizer::getGSLMassMatrix(
std::vector<complex> const &massMatrix
) const
{
gsl_matrix_complex *matrix = gsl_matrix_complex_alloc(m_N, m_N);
for (size_t i = 0; i != m_N; ++i)
for (size_t j = 0; j != m_N; ++j)
gsl_matrix_complex_set(
matrix,
i, j,
gsl_complex_rect(REAL(massMatrix[index(i, j)]),
IMAG(massMatrix[index(i, j)]))
);
return matrix;
}
void Diagonalizer::loadGSLMatrix(
gsl_matrix_complex const *matrix,
std::vector<complex> &target
)
{
target = std::vector<complex>(m_N * m_N);
for (size_t i = 0; i != m_N; ++i)
for (size_t j = 0; j != m_N; ++j) {
auto gslcomplex = gsl_matrix_complex_get(matrix, i, j);
target[index(i, j)] = {GSL_REAL(gslcomplex),
GSL_IMAG(gslcomplex)};
}
}
void Diagonalizer::loadGSLVector(
gsl_vector const *vect,
std::vector<real> &target
)
{
target = std::vector<real>(m_N);
for (size_t i = 0; i != m_N; ++i)
target[i] = gsl_vector_get(vect, i);
}
void Diagonalizer::positiveDiagonal(
std::vector<complex> &transfer
)
{
for (size_t j = 0; j != m_N; ++j) {
if (REAL(transfer[index(j, j)]) < 0) {
for (size_t i = 0; i != m_N; ++i)
transfer[index(i, j)] *= -1;
}
}
}
std::vector<Diagonalizer::complex> Diagonalizer::hermitian(
std::vector<complex> const &init
) const
{
std::vector<complex> res(init);
for (size_t i = 0; i != m_N; ++i)
for (size_t j = 0; j != m_N; ++j)
res[index(i, j)] = CONJ(init[index(j, i)]);
return res;
}
std::vector<Diagonalizer::complex> Diagonalizer::dot(
std::vector<complex> const &A,
std::vector<complex> const &B
) const
{
std::vector<complex> C(A);
for (size_t i = 0; i != m_N; ++i)
for (size_t j = 0; j != m_N; ++j) {
C[index(i, j)] = 0;
for (size_t k = 0; k != m_N; ++k)
C[index(i, j)] += A[index(i, k)] * B[index(k, j)];
}
return C;
}
......@@ -886,8 +886,8 @@ namespace csl {
vec.erase(last, vec.end());
};
for (auto &diag : diagData) {
cutSimilar(diag.mixings);
cutSimilar(diag.masses);
// cutSimilar(diag.mixings);
// cutSimilar(diag.masses);
cutSimilar(diag.dependencies);
}
// for (auto &f : functions) {
......@@ -1173,9 +1173,9 @@ namespace csl {
diag.masses.begin(),
diag.masses.end()
);
std::sort(res.begin(), res.end());
auto last = std::unique(res.begin(), res.end());
res.erase(last, res.end());
// std::sort(res.begin(), res.end());
// auto last = std::unique(res.begin(), res.end());
// res.erase(last, res.end());
return res;
}
......
......@@ -114,7 +114,7 @@ namespace mty::mssm_input {
/**
* @brief Tau trilinear coupling.
*/
inline csl::Expr Atau = csl::constant_s("A_\\tau");
inline csl::Expr Atau = csl::constant_s("A_tau ; A_{\\tau}");
/**
* @brief Down-type Higgs mass squared.
......
......@@ -279,7 +279,7 @@ class GaugedGroup {
bool broken;
csl::Expr coupling;
csl::Expr coupling;
std::string nameGenerator = "T";
......
......@@ -219,7 +219,8 @@ public:
WilsonSet computeWilsonCoefficients(
int order,
std::vector<Insertion> const &insertions,
FeynOptions feynOptions = {}
FeynOptions feynOptions = {},
bool disableFermionOrdering = false
);
WilsonSet computeWilsonCoefficients_default(
......
......@@ -354,6 +354,9 @@ public:
std::vector<mty::Lagrangian::TermType> &interaction
);
void gatherMass(mty::Particle const &part);
void gatherMass(std::string const &name);
void gatherMasses();
void refresh();
......
......@@ -241,6 +241,7 @@ protected:
csl::Expr V_us;
csl::Expr A_b;
csl::Expr A_t;
csl::Expr A_tau;
csl::Expr U_sb_00;
csl::Expr N_u4;
csl::Expr V_ts;
......
......@@ -211,12 +211,6 @@ protected:
void initInteractions162();
void initInteractions163();
void initInteractions164();
void initInteractions165();
void initInteractions166();
void initInteractions167();
void initInteractions168();
void initInteractions169();
void initInteractions170();
void initSpectrum();
protected:
......
......@@ -660,6 +660,7 @@ void PMSSM_LEM::initContent()
A_t = mssm_input::At;
A_b = mssm_input::Ab;
A_tau = mssm_input::Atau;
e_em = sm_input::e_em;
mu_h = mssm_input::mu;
V_ub_mod = sm_input::V_ub_mod;
......@@ -921,12 +922,6 @@ void PMSSM_LEM::initInteractions()
initInteractions162();
initInteractions163();
initInteractions164();
initInteractions165();
initInteractions166();
initInteractions167();
initInteractions168();
initInteractions169();
initInteractions170();
}
void PMSSM_LEM::initSpectrum()
......@@ -1014,7 +1009,7 @@ void PMSSM_LEM::initSpectrum()
, csl::pow_s(csl::sin_s(beta), 2)})})}) , csl::prod_s({csl::intfraction_s(1, 2), csl::pow_s(M_Z, 2), csl::pow_s(csl::sin_s(beta), 2), csl::pow_s(csl::sin_s(theta_W), 2)}) , csl::prod_s({csl::intfraction_s(-1, 2), csl::pow_s(M_Z, 2), csl::sum_s({1 , csl::prod_s({-1
, csl::pow_s(csl::sin_s(beta), 2)})}), csl::pow_s(csl::sin_s(theta_W), 2)})}), csl::sum_s({csl::prod_s({-2
, mu_h, m_tau, csl::pow_s(csl::cos_s(beta), (-1)), csl::sin_s(beta)}) , csl::prod_s({2
, csl::pow_s(2, csl::intfraction_s(1, 2)), A_b, M_W, csl::pow_s(e_em, (-1)), csl::cos_s(beta), csl::sin_s(theta_W)})})},
, csl::pow_s(2, csl::intfraction_s(1, 2)), A_tau, M_W, csl::pow_s(e_em, (-1)), csl::cos_s(beta), csl::sin_s(theta_W)})})},
{0, csl::sum_s({csl::pow_s(M_tauR, 2) , csl::prod_s({csl::pow_s(m_tau, 2), csl::pow_s(csl::cos_s(beta), (-2)), csl::sum_s({1 , csl::prod_s({-1
, csl::pow_s(csl::sin_s(beta), 2)})})}) , csl::prod_s({csl::pow_s(M_Z, 2), csl::pow_s(csl::sin_s(beta), 2), csl::pow_s(csl::sin_s(theta_W), 2)}) , csl::prod_s({-1
, csl::pow_s(M_Z, 2), csl::sum_s({1 , csl::prod_s({-1
......
This diff is collapsed.
......@@ -75,8 +75,7 @@ namespace sgl {
// If double application, we continue only if the two
// chains are composed of the fermions they want AND are
// canonical
if (twice
&& areCanonical(c1, c2)
if (areCanonical(c1, c2)
&& c1.isHappy()
&& c2.isHappy())
continue;
......@@ -88,8 +87,10 @@ namespace sgl {
GExpr res;
if (twice)
res = c1.applyGeneralFierzTwice(c2);
else
else if (!c1.isHappy() || !c2.isHappy())
res = c1.applyGeneralFierz(c2);
else
continue;
sub[i] = res;
sub[j] = cslexpr_s(1);
return true;
......
......@@ -138,7 +138,7 @@ void PMSSM_Model::approximateInputMatrices()
Au->setTensor(A1);
Ad->setTensor(A2);
Ae->setTensor(A2);
Ae->setTensor(A3);
}
void PMSSM_Model::approximateQuarkMasses()
......
......@@ -30,7 +30,8 @@ SM_Model::SM_Model(bool initialize)
:mty::Model("models/files/SM.json")
{
getParticle("G")->setDrawType(drawer::ParticleType::Gluon);
replace(getScalarCoupling("g_s"), sm_input::g_s);
replace(getScalarCoupling("g_s"), g_s);
getGaugedGroup("SU3c")->setCouplingConstant(g_s);
Particle H = GetParticle(*this, "H");
......@@ -137,9 +138,11 @@ void SM_Model::diagonalizeSMMassMatrices()
// Diagonalizing what can be
///////////////////////////////////////////////////
DiagonalizeMassMatrices(*this);
diagonalizeMassMatrices();
renameParticle("B", "A");
renameParticle("W_3", "Z");
gatherMass("W");
gatherMass("Z");
csl::Expr gY = getScalarCoupling("g_Y");
csl::Expr gL = getScalarCoupling("g_L");
......
......@@ -534,10 +534,13 @@ WilsonSet Model::getWilsonCoefficients(
WilsonSet Model::computeWilsonCoefficients(
int order,
std::vector<Insertion> const &insertions,
FeynOptions feynOptions
FeynOptions feynOptions,
bool disableFermionOrdering
)
{
feynOptions.orderExternalFermions = true;
if (!disableFermionOrdering) {
feynOptions.orderExternalFermions = true;
}
if (order == TreeLevel) {
return computeWilsonCoefficients_default(
order, insertions, feynOptions
......
......@@ -834,14 +834,49 @@ void ModelBuilder::applyUnitaryCondition(
std::vector<std::vector<csl::Expr>> const &unitary
)
{
for (size_t i = 0; i != unitary.size(); ++i)
for (size_t j = 0; j != unitary.size(); ++j) {
csl::Expr expr = CSL_0;
for (size_t k = 0; k != unitary.size(); ++k)
expr += csl::GetComplexConjugate(unitary[k][i])*unitary[k][j];
replace(expr, (i == j) ? CSL_1 : CSL_0);
// replace(csl::GetComplexConjugate(expr), (i == j) ? CSL_1 : CSL_0);
auto dependsOnMixing = [&](csl::Expr const &expr) {
for (const auto &row : unitary)
for (const auto el : row)
if (expr->dependsExplicitlyOn(el.get()))
return true;
return false;
};
std::vector<csl::Expr> terms = clearDependencies(
L.interaction,
[&](Lagrangian::TermType const &term) {
return dependsOnMixing(term->getTerm());
});
for (auto &term : terms) {
bool applied = false;
for (size_t i = 0; i != unitary.size(); ++i) {
for (size_t j = 0; j != unitary.size(); ++j) {
csl::Expr lhs = GetComplexConjugate(unitary[0][i]);
csl::Expr rhs = CSL_0;
for (size_t k = 1; k != unitary.size(); ++k)
rhs += csl::GetComplexConjugate(unitary[k][i])*unitary[k][j];
csl::Expr remnant = (i == j) ? 1 : 0;
rhs = (remnant - rhs) / unitary[0][j];
csl::Expr replaced = csl::Replaced(term, lhs, rhs);
csl::DeepExpand(replaced);
csl::DeepHardFactor(replaced);
if (!dependsOnMixing(replaced)) {
applied = true;
term = replaced;
break;
}
}
if (applied) break;
}