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

Right-handed charges corrected

parent fcedd5d9
......@@ -29,9 +29,9 @@ marty:
install:
cd marty&& $(MAKE) install
cd csl&& $(MAKE) install
cd grafed&& $(MAKE) install
cd jsonparser&& $(MAKE) install
cd dependencies&& $(MAKE) install
cd grafed&& $(MAKE) install
@echo "$(GREEN)$(BOLD)The installation is finished, you may now watch the trilogy 'Back to the future' "
@echo "to honor the McFly family !"
@echo "(or start using the code, as you want)$(NC)$(NORMAL)"
......@@ -39,6 +39,6 @@ install:
clean:
cd marty&& $(MAKE) clean;
cd csl&& $(MAKE) clean;
cd grafed&& $(MAKE) clean;
cd jsonparser&& $(MAKE) clean;
cd dependencies&& $(MAKE) clean;
cd grafed&& $(MAKE) clean;
......@@ -74,6 +74,7 @@
#include "timeMonitor.h"
#include "linear_map.h"
#include "initSanitizer.h"
#include "multipermutation.h"
#include "librarygenerator.h"
#include "libraryfunction.h"
......
......@@ -4,7 +4,7 @@ namespace csl {
class Expr;
int matchBOnA(csl::Expr const& A, csl::Expr &B);
int matchBOnA(csl::Expr const& A, csl::Expr const &B);
bool hardComparison(csl::Expr const&, csl::Expr const&);
bool hardOrdering(csl::Expr const&, csl::Expr const&);
......
#pragma once
#include <vector>
namespace csl {
class MultiPermutation {
public:
MultiPermutation(std::vector<std::vector<size_t>> const &t_numbers);
MultiPermutation(MultiPermutation const &other) = delete;
MultiPermutation &operator=(MultiPermutation const &other) = delete;
MultiPermutation(MultiPermutation &&other) = default;
MultiPermutation &operator=(MultiPermutation &&other) = default;
bool nextPermutation();
template<class T>
std::vector<T> applyPermutation(std::vector<T> const &collection)
{
std::vector<T> res(collection);
for (size_t i = 0; i != m_numbers.size(); ++i) {
assert(m_numbers_init[i] < collection.size()
&& m_numbers[i] < collection.size());
res[m_numbers_init[i]] = collection[m_numbers[i]];
}
return res;
}
void print() const;
private:
void load(std::vector<std::vector<size_t>> const &t_numbers);
private:
std::vector<size_t> m_numbers;
std::vector<size_t> m_numbers_init;
std::vector<std::vector<size_t>::iterator> m_boundaries;
};
}
......@@ -401,17 +401,20 @@ std::optional<Expr> Abbrev::findExisting(
if (pos < abbreviations.size()
&& abbreviations[pos]->getEncapsulated()->compareWithDummy(
encapsulated.get())) {
return dynamic_cast<LiteralParent*>(abbreviations[pos])->generateInstance();
return dynamic_cast<LiteralParent*>(abbreviations[pos])
->generateInstance();
}
if (pos-1 < abbreviations.size()
&& abbreviations[pos-1]->getEncapsulated()->compareWithDummy(
encapsulated.get())) {
return dynamic_cast<LiteralParent*>(abbreviations[pos-1])->generateInstance();
return dynamic_cast<LiteralParent*>(abbreviations[pos-1])
->generateInstance();
}
if (pos+1 < abbreviations.size()
&& abbreviations[pos+1]->getEncapsulated()->compareWithDummy(
encapsulated.get())) {
return dynamic_cast<LiteralParent*>(abbreviations[pos+1])->generateInstance();
return dynamic_cast<LiteralParent*>(abbreviations[pos+1])
->generateInstance();
}
}
else {
......@@ -461,6 +464,8 @@ std::optional<Expr> Abbrev::findExisting(
std::map<csl::Index, csl::Index> mapping;
if (encapsulated->compareWithDummy(comparison.get(), mapping))
return (*ab)(structure.getIndex());
if (csl::hardComparison(encapsulated, comparison))
return (*ab)(structure.getIndex());
}
}
......
......@@ -262,7 +262,7 @@ void VisitEachNodeCut(
return;
if (depth != 0) {
for (size_t i = 0; i != init->size(); ++i)
VisitEachNode((*init)[i], f, depth-1);
VisitEachNodeCut((*init)[i], f, depth-1);
}
}
......@@ -276,7 +276,7 @@ void VisitEachNodeCut(
return;
if (depth != 0) {
for (size_t i = 0; i != init->size(); ++i)
VisitEachNode((*init)[i].get(), f, depth-1);
VisitEachNodeCut((*init)[i].get(), f, depth-1);
}
}
......
#include <numeric>
#include "interface.h"
#include "abreviation.h"
#include "multipermutation.h"
namespace csl {
static MultiPermutation getAllPossiblePermutations(
std::vector<csl::Expr> const &tensors
)
{
std::vector<std::vector<size_t>> permutations;
permutations.reserve(tensors.size());
std::vector<size_t> indicesLeft(tensors.size());
std::iota(indicesLeft.begin(), indicesLeft.end(), 0);
while (!indicesLeft.empty()) {
auto const &tensor1 = tensors[indicesLeft[0]];
permutations.push_back({indicesLeft[0]});
indicesLeft.erase(indicesLeft.begin());
if (!tensor1->getCommutable())
continue;
for (size_t i = 0; i != indicesLeft.size(); ++i) {
auto const &tensor2 = tensors[indicesLeft[i]];
if (tensor1->getParent_info() == tensor2->getParent_info()
&& tensor1->isComplexConjugate() == tensor2->isComplexConjugate()) {
permutations.back().push_back(i);
indicesLeft.erase(indicesLeft.begin() + i);
--i;
}
}
}
return MultiPermutation { permutations };
}
static void sortTensors(std::vector<csl::Expr> &tensors)
{
auto free = [&](csl::IndexStructure const &index) {
......@@ -42,25 +71,41 @@ static void sortTensors(std::vector<csl::Expr> &tensors)
tensors = std::move(sorted);
}
int matchBOnA(csl::Expr const& A, csl::Expr &B)
static std::pair<std::vector<csl::Expr>, std::vector<csl::Expr>>
getSortedTensors(csl::Expr const &A, csl::Expr const &B)
{
std::vector<csl::Expr> tensorsInA;
std::vector<csl::Expr> tensorsInB;
csl::VisitEachLeaf(A, [&](csl::Expr const& el)
csl::VisitEachNodeCut(A, [&](csl::Expr const& el)
{
if (el->isIndexed())
if (csl::IsIndicialTensor(el))
tensorsInA.push_back(el);
return csl::IsSum(el);
});
csl::VisitEachLeaf(B, [&](csl::Expr const& el)
csl::VisitEachNodeCut(B, [&](csl::Expr const& el)
{
if (el->isIndexed())
if (csl::IsIndicialTensor(el))
tensorsInB.push_back(el);
return csl::IsSum(el);
});
// if (tensorsInA.size() != tensorsInB.size()) {
// return tensorsInA.size() < tensorsInB.size();
// }
sortTensors(tensorsInA);
sortTensors(tensorsInB);
return {tensorsInA, tensorsInB};
}
int matchBOnA(
csl::Expr &B,
std::vector<csl::Expr> const &tensorsInA,
std::vector<csl::Expr> const &tensorsInB
)
{
if (tensorsInA.size() != tensorsInB.size()) {
return tensorsInA.size() < tensorsInB.size();
}
sortTensors(tensorsInA);
sortTensors(tensorsInB);
std::vector<std::pair<csl::Index, csl::Index>> mapping;
for (size_t i = 0; i != tensorsInA.size(); ++i) {
if (tensorsInA[i]->getParent_info()
......@@ -127,15 +172,35 @@ int matchBOnA(csl::Expr const& A, csl::Expr &B)
return -1;
}
int matchBOnA(csl::Expr const& A, csl::Expr &B)
{
auto [tensorsInA, tensorsInB] = getSortedTensors(A, B);
if (tensorsInA.size() != tensorsInB.size()) {
return tensorsInA.size() < tensorsInB.size();
}
return matchBOnA(B, tensorsInA, tensorsInB);
}
static bool hardComparison_impl(
csl::Expr const &A,
csl::Expr &B)
csl::Expr const &B)
{
const int match = matchBOnA(A, B);
if (match != -1)
return false;
auto res = A->compareWithDummy(B.get());
return res;
auto [tensorsInA, tensorsInB] = getSortedTensors(A, B);
MultiPermutation permutation = getAllPossiblePermutations(tensorsInB);
size_t i = 0;
do {
csl::Expr B_cpy = csl::DeepCopy(B);
auto tensorsInB_perm = permutation.applyPermutation(tensorsInB);
const int match = matchBOnA(B_cpy, tensorsInA, tensorsInB_perm);
if (match != -1)
return false;
if (A->compareWithDummy(B_cpy.get()))
return true;
++i;
} while (permutation.nextPermutation());
return false;
}
bool hardComparison(
......
#include "multipermutation.h"
#include <algorithm>
#include <iostream>
namespace csl {
MultiPermutation::MultiPermutation(
std::vector<std::vector<size_t>> const &t_numbers
)
{
load(t_numbers);
}
void MultiPermutation::load(
std::vector<std::vector<size_t>> const &t_numbers
)
{
size_t sz = 0;
for (const auto &vec : t_numbers)
sz += vec.size();
m_numbers.reserve(sz);
m_boundaries.reserve(t_numbers.size());
for (const auto &vec : t_numbers) {
m_boundaries.push_back(m_numbers.end());
m_numbers.insert(m_numbers.end(), vec.begin(), vec.end());
}
m_boundaries.push_back(m_numbers.end());
m_numbers_init = m_numbers;
}
bool MultiPermutation::nextPermutation()
{
for (size_t i = 0; i+1 != m_boundaries.size(); ++i) {
if (std::next_permutation(m_boundaries[i], m_boundaries[i+1]))
return true;
}
return false;
}
void MultiPermutation::print() const
{
for (size_t i = 0; i+1 != m_boundaries.size(); ++i) {
std::cout << "( ";
for (auto iter = m_boundaries[i]; iter != m_boundaries[i+1];
++iter) {
std::cout << *iter << " ";
}
std::cout << ")" << '\n';
}
}
}
......@@ -80,7 +80,7 @@ namespace mty::mssm_input {
/**
* @brief mu parameter.
*/
inline csl::Expr mu = csl::constant_s("mu");
inline csl::Expr mu = csl::constant_s("mu_h");
/**
* @brief Sign of the bilinear Higgs term in the superpotential.
*/
......@@ -212,57 +212,57 @@ namespace mty::mssm_input {
*/
inline csl::Expr MqbR = csl::constant_s("M_qbR");
#define CST(arg) csl::constant_s(arg)
#define MAT(name) \
#define MTY_CST(arg) csl::constant_s(arg)
#define MTY_MAT(name) \
csl::matrix_s({\
{CST(#name"_11"), CST(#name"_12"), CST(#name"_13")},\
{CST(#name"_21"), CST(#name"_22"), CST(#name"_23")},\
{CST(#name"_31"), CST(#name"_32"), CST(#name"_33")}\
{MTY_CST(#name"_11"), MTY_CST(#name"_12"), MTY_CST(#name"_13")},\
{MTY_CST(#name"_21"), MTY_CST(#name"_22"), MTY_CST(#name"_23")},\
{MTY_CST(#name"_31"), MTY_CST(#name"_32"), MTY_CST(#name"_33")}\
})
#define MAT_SPEC(name, u, c, t) \
#define MTY_MAT_SPEC(name, u, c, t) \
csl::matrix_s({\
{CST(#name"_"#u#u), CST(#name"_"#u#c), CST(#name"_"#u#t)},\
{CST(#name"_"#c#u), CST(#name"_"#c#c), CST(#name"_"#c#t)},\
{CST(#name"_"#t#u), CST(#name"_"#t#c), CST(#name"_"#t#t)}\
{MTY_CST(#name"_"#u#u), MTY_CST(#name"_"#u#c), MTY_CST(#name"_"#u#t)},\
{MTY_CST(#name"_"#c#u), MTY_CST(#name"_"#c#c), MTY_CST(#name"_"#c#t)},\
{MTY_CST(#name"_"#t#u), MTY_CST(#name"_"#t#c), MTY_CST(#name"_"#t#t)}\
})
/**
* @brief Left squark soft SUSY-brekaing mass matrix.
*/
inline csl::Expr MSQ2 = MAT(MSQ2);
inline csl::Expr MSQ2 = MTY_MAT(MSQ2);
/**
* @brief Left slepton soft SUSY-brekaing mass matrix.
*/
inline csl::Expr MSL2 = MAT(MSL2);
inline csl::Expr MSL2 = MTY_MAT(MSL2);
/**
* @brief Right up-type squark soft SUSY-brekaing mass matrix.
*/
inline csl::Expr MSu2 = MAT_SPEC(MSu2, u, c, t);
inline csl::Expr MSu2 = MTY_MAT_SPEC(MSu2, u, c, t);
/**
* @brief Right down-type squark soft SUSY-brekaing mass matrix.
*/
inline csl::Expr MSd2 = MAT_SPEC(MSd2, d, s, b);
inline csl::Expr MSd2 = MTY_MAT_SPEC(MSd2, d, s, b);
/**
* @brief Right selectron soft SUSY-brekaing mass matrix.
*/
inline csl::Expr MSe2 = MAT_SPEC(MSe2, e, mu, tau);
inline csl::Expr MSe2 = MTY_MAT_SPEC(MSe2, e, mu, tau);
/**
* @brief Up-type squark soft SUSY-breaking trilinear couplings
*/
inline csl::Expr Tu = MAT_SPEC(T_u, u, c, t);
inline csl::Expr Tu = MTY_MAT_SPEC(T_u, u, c, t);
/**
* @brief Down-type squark soft SUSY-breaking trilinear couplings
*/
inline csl::Expr Td = MAT_SPEC(T_d, d, s, b);
inline csl::Expr Td = MTY_MAT_SPEC(T_d, d, s, b);
/**
* @brief Slepton soft SUSY-breaking trilinear couplings
*/
inline csl::Expr Te = MAT_SPEC(T_e, e, mu, tau);
inline csl::Expr Te = MTY_MAT_SPEC(T_e, e, mu, tau);
#undef CST
#undef MAT
#undef MAT_SPEC
#undef MTY_CST
#undef MTY_MAT
#undef MTY_MAT_SPEC
}
namespace mty {
......
......@@ -209,8 +209,6 @@ protected:
void initInteractions160();
void initInteractions161();
void initInteractions162();
void initInteractions163();
void initInteractions164();
void initSpectrum();
protected:
......
......@@ -920,8 +920,6 @@ void PMSSM_LEM::initInteractions()
initInteractions160();
initInteractions161();
initInteractions162();
initInteractions163();
initInteractions164();
}
void PMSSM_LEM::initSpectrum()
......
This diff is collapsed.
......@@ -404,7 +404,7 @@ void MSSM_Model::initGauginoInteractions(
csl::Expr scal = (isLeft) ? cc(scalar(s_index)) : scalar(s_index);
csl::Expr ferm = (isLeft) ? fermion(f_index) : cc(fermion(f_index));
csl::Expr connexion = (isLeft) ?
dirac4.C_matrix({a, b}) : -dirac4.getDelta()({a, b});
dirac4.C_matrix({a, b}) : dirac4.getDelta()({a, b});
addLagrangianTerm(
csl::DeepRefreshed(csl::prod_s({
-csl::sqrt_s(2), coupling, T({A, i, j}),
......@@ -441,6 +441,8 @@ void MSSM_Model::initU1GauginoInteractions(
bool isLeft = (fermion->getChirality() == Chirality::Left);
csl::Expr scal = (isLeft) ? cc(scalar(s_index)) : scalar(s_index);
csl::Expr ferm = (isLeft) ? fermion(f_index) : cc(fermion(f_index));
// Minus sign here if right-handed because the charge changes sign,
// this sign is not present for non-abelian couplings T^A
csl::Expr connexion = (isLeft) ?
dirac4.C_matrix({a, b}) : -dirac4.getDelta()({a, b});
addLagrangianTerm(
......@@ -640,29 +642,29 @@ void MSSM_Model::initTrilinears()
// Up-type tri-linear from the super potential
addLagrangianTerm(
mu_h
* cc(s_Qi({I, A, i}))
* cc(Yu({J, I}))
* s_Ui({J, A})
* Hd(i),
cc(mu_h)
* cc(s_Ui({J, A}))
* Yu({J, I})
* s_Qi({I, A, i})
* cc(Hd(i)),
true // Add also the complex conjugate of this term
);
// Down-type tri-linear from the super potential
addLagrangianTerm(
mu_h
* cc(s_Qi({I, A, i}))
* cc(Yd({J, I}))
* s_Di({J, A})
* Hu(i),
cc(mu_h)
* cc(s_Di({J, A}))
* Yd({J, I})
* s_Qi({I, A, i})
* cc(Hu(i)),
true // Add also the complex conjugate of this term
);
// Lepton tri-linear from the super potential
addLagrangianTerm(
mu_h
* cc(s_Li({I, i}))
* cc(Ye({J, I}))
* s_Ei({J})
* Hu(i),
cc(mu_h)
* cc(s_Ei({J}))
* Ye({J, I})
* s_Li({I, i})
* cc(Hu(i)),
true // Add also the complex conjugate of this term
);
}
......@@ -681,6 +683,8 @@ void MSSM_Model::initQuarticDTerms()
== mty::ParticleType::GhostBoson)
continue;
auto rep1 = particles[p1]->getGroupIrrep((*gauge)[i]);
bool sign1 = SUSY[particles[p1]]->isFermionic()
&& SUSY[particles[p1]]->getChirality() == Chirality::Right;
for (size_t p2 = 0; p2 < particles.size(); ++p2) {
if (particles[p2]->getSpinDimension() != 1
or particles[p2]->getParticleType()
......@@ -690,6 +694,8 @@ void MSSM_Model::initQuarticDTerms()
if (rep1 != rep2
and (*gauge)[i]->getType() != mty::group::Type::U1)
continue;
bool sign2 = SUSY[particles[p2]]->isFermionic()
&& SUSY[particles[p2]]->getChirality() == Chirality::Right;
csl::Index mu = mty::MinkowskiIndex();
csl::Expr field1 = particles[p1]->getInstance();
......@@ -729,8 +735,9 @@ void MSSM_Model::initQuarticDTerms()
}
});
csl::Replace(quadratic2, B, A);
csl::Expr sign = (sign1 ^ sign2) ? CSL_M_1 : CSL_1;
addLagrangianTerm(
csl::Refreshed(CSL_HALF
csl::Refreshed(sign * CSL_HALF
* quadratic1 * quadratic2
));
}
......@@ -1097,10 +1104,11 @@ void MSSM_Model::diagonalize2By2Matrices()
csl::Expr mu_2 = cc(mu_h) * mu_h;
csl::Expr MZ2 = csl::pow_s(sm_input::M_Z, 2);
Replaced(*this,
g_L*g_L + g_Y*g_Y,
4 * MZ2 / (v_h*v_h));
csl::Expr c2beta = 1 - 2 * csl::pow_s(csl::sin_s(beta_h), 2);
// Replaced(*this,
// g_L*g_L + g_Y*g_Y,
// csl::pow_s(sm_input::e_em / csl::sin_s(sm_input::theta_W), 2)
// * (1 + csl::pow_s(csl::tan_s(sm_input::theta_W), 2)));
csl::Expr c2beta = csl::cos_s(2*beta_h);
Replaced(*this,
m_sHu,
csl::sqrt_s(M_A0*M_A0 - m_sHd*m_sHd - 2*mu_2));
......@@ -1110,10 +1118,10 @@ void MSSM_Model::diagonalize2By2Matrices()
csl::Expanded(-MZ2 * c2beta / 2)
+ csl::Expanded(M_A0*M_A0 * (1 - c2beta) / 2)
- mu_2));
Replaced(*this,
csl::pow_s(csl::cos_s(beta_h), 2),
1 - csl::pow_s(csl::sin_s(beta_h), 2)
);
// Replaced(*this,
// csl::pow_s(csl::cos_s(beta_h), 2),
// 1 - csl::pow_s(csl::sin_s(beta_h), 2)
// );
mty::Particle eta_u = getParticle("eta_u");
mty::Particle eta_d = getParticle("eta_d");
mty::Particle G0 = scalarboson_s("G0; G^0", *this);
......@@ -1263,9 +1271,9 @@ void MSSM_Model::diagonalizeYukawas()
}
void MSSM_Model::adjustCouplingConstants()
{
csl::Expr c2 = csl::pow_s(csl::cos_s(beta_h), 2);
csl::Expr s2 = csl::pow_s(csl::sin_s(beta_h), 2);
replace(c2, 1 - s2);
//csl::Expr c2 = csl::pow_s(csl::cos_s(beta_h), 2);
//csl::Expr s2 = csl::pow_s(csl::sin_s(beta_h), 2);
//replace(c2, 1 - s2);
csl::Expr e = sm_input::e_em;
csl::Expr thetaW = sm_input::theta_W;
......@@ -1278,9 +1286,9 @@ void MSSM_Model::adjustCouplingConstants()
replace(
g_L,
e / csl::sin_s(thetaW));
replace(
csl::cos_s(thetaW),
sm_input::M_W / sm_input::M_Z);
// replace(
// csl::cos_s(thetaW),
// sm_input::M_W / sm_input::M_Z);
}
void MSSM_Model::breakSMFlavorSymmetry()
{
......@@ -1541,8 +1549,8 @@ MSSM_HEM::MSSM_HEM(
promoteToMajorana("sG");
std::cout << "Gathering MSSM inputs ..." << std::endl;
gatherMSSMInputs();
std::cout << "Checking Hermiticity ..." << std::endl;
checkHermiticity();
// std::cout << "Checking Hermiticity ..." << std::endl;
// checkHermiticity();
computeFeynmanRules();
if (save) {
......
......@@ -1062,8 +1062,11 @@ auto replaceMajorana(
const size_t sz = expr->size();
bool left = false;
auto d = dirac4.getDelta();
auto P_L = dirac4.P_L;
auto P_R = dirac4.P_R;
auto P_1 = dirac4.P_L;
auto P_2 = dirac4.P_R;
if (xi->getChirality() == Chirality::Right) {
std::swap(P_1, P_2);
}
auto cc = csl::GetComplexConjugate;
auto C = dirac4.C_matrix;
bool leftMajorana = false;
......@@ -1085,26 +1088,26 @@ auto replaceMajorana(
const csl::Expr f = (leftMajorana) ? CSL_HALF : CSL_1;
if (left && conjug) {
const auto P = (alreadyProjected) ?
f*d({b, a}) : P_R({b, a});
f*d({b, a}) : P_2({b, a});