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

Signs of conjugation && Fermion FR corrected

parent 1a46b7cc
......@@ -4228,6 +4228,27 @@ std::vector<Expr> sumDummyIndices(
return res;
}
void sortIndices(
std::vector<csl::Index> &indices,
std::vector<size_t> &positions
)
{
for (size_t i = 0; i != indices.size(); ++i) {
size_t posMini = i;
for (size_t j = i+1; j != indices.size(); ++j) {
if (indices[j] < indices[posMini])
posMini = j;
}
if (posMini != i) {
std::swap(indices[i], indices[posMini]);
std::for_each(positions.begin(), positions.end(), [&](size_t &s) {
if (s == posMini) s = i;
else if (s == i) s = posMini;
});
}
}
}
csl::vector_expr IProd::breakSpace(
const Space* brokenSpace,
const vector<const Space*>& newSpaces,
......@@ -4265,6 +4286,7 @@ csl::vector_expr IProd::breakSpace(
std::distance(brokenIndices.begin(), pos)
);
}
sortIndices(brokenIndices, dummyIndices);
vector<csl::vector_expr> brokenExpr(argument.size());
for (size_t i = 0; i != argument.size(); ++i) {
......
......@@ -10,7 +10,7 @@ DYNAMICLIBRARY = libmarty.so
LTVERSION = 2.15
CXXDEBUG = $(CXX) -g
CXXFLAGS = -Wall -Wextra -Wpedantic -std=c++17 -O3 -fPIC $(MARTY_CXX_FLAGS)
CXXFLAGS = -Wall -Wextra -Wpedantic -std=c++17 -O3 -fPIC $(MARTY_CXX_FLAGS)
CXXNOOPTIFLAGS = -Wall -Wextra -Wpedantic -std=c++17 -fPIC $(MARTY_CXX_FLAGS)
SRCDIR = src
......
......@@ -357,6 +357,15 @@ private:
std::string saveFile;
};
class MSSM_HEM: public MSSM_Model {
public:
MSSM_HEM(
std::string const &flhaFile = "",
std::string const &saveFile = ""
);
};
} // End of namespace mty
#endif
#pragma once
#include "PMSSM.h"
namespace mty {
class NMFV_Model: public PMSSM_Model {
public:
enum class Mixing {
QL, UR, DR, LL, ER, Au, Ad, Ae
};
struct MixingInfo {
Mixing mixing;
size_t i;
size_t j;
};
NMFV_Model(
std::vector<MixingInfo> const &mixings,
std::string const &slhaFile = "",
std::string const &saveFile = ""
);
NMFV_Model(
std::initializer_list<MixingInfo> mixings,
std::string const &slhaFile = "",
std::string const &saveFile = ""
)
:NMFV_Model(
std::vector<MixingInfo>(mixings),
slhaFile,
saveFile
)
{
}
void addMixing(MixingInfo mix);
};
std::ostream &operator<<(
std::ostream &out,
NMFV_Model::Mixing mixing
);
}
......@@ -33,11 +33,10 @@ public:
PMSSM_Model(
std::string const &slhaFile = "",
std::string const &saveFile = ""
std::string const &saveFile = "",
bool init = true
);
~PMSSM_Model();
protected:
void approximateYukawa();
......@@ -47,6 +46,14 @@ protected:
void approximateSFermionMixings();
void renameSFermions();
void getToLowEnergyLagrangian();
void addAllowedMixing(std::vector<std::string> const &names);
void mergeAllowedMixings();
bool isSuppressedMixing(mty::InteractionTerm const &massTerm) const;
protected:
std::vector<std::set<std::string>> allowedMixings;
};
} // End of namespace mty
......
......@@ -204,61 +204,61 @@ namespace mty::sm_input {
);
///////////////////////////////////////////////////
// CKM Matrix
// CKM Matrix PDG 2020
///////////////////////////////////////////////////
inline
csl::Expr V_ud_mod = csl::constant_s(
"V_ud_mod",
csl::float_s(0.97420)
csl::float_s(0.97370)
);
inline
csl::Expr V_us_mod = csl::constant_s(
"V_us_mod",
csl::float_s(0.2243)
csl::float_s(0.2245)
);
inline
csl::Expr V_ub_mod = csl::constant_s(
"V_ub_mod",
csl::float_s(3.94e-03)
csl::float_s(3.82e-03)
);
inline
csl::Expr V_cb_mod = csl::constant_s(
"V_cb_mod",
csl::float_s(42.2e-03)
csl::float_s(41.0e-03)
);
inline
csl::Expr V_cd_mod = csl::constant_s(
"V_cd_mod",
csl::float_s(0.2180)
csl::float_s(0.2210)
);
inline
csl::Expr V_cs_mod = csl::constant_s(
"V_cs_mod",
csl::float_s(0.997)
csl::float_s(0.987)
);
inline
csl::Expr V_td_mod = csl::constant_s(
"V_td_mod",
csl::float_s(8.1e-03)
csl::float_s(8.0e-03)
);
inline
csl::Expr V_ts_mod = csl::constant_s(
"V_ts_mod",
csl::float_s(39.4e-03)
csl::float_s(38.8e-03)
);
inline
csl::Expr V_tb_mod = csl::constant_s(
"V_tb_mod",
csl::float_s(1.019)
csl::float_s(1.013)
);
inline
......@@ -296,29 +296,29 @@ namespace mty::sm_input {
inline
csl::Expr V_ud = csl::constant_s(
"V_ud",
csl::float_s(0.97420),
V_ud_mod->getValue(),
csl::ComplexProperty::Real
);
inline
csl::Expr V_us = csl::constant_s(
"V_us",
csl::float_s(0.2243),
V_us_mod->getValue(),
csl::ComplexProperty::Real
);
inline
csl::Expr V_cb = csl::constant_s(
"V_cb",
csl::float_s(42.2e-03),
V_cb_mod->getValue(),
csl::ComplexProperty::Real
);
inline
csl::Expr V_tb = csl::constant_s(
"V_tb",
csl::float_s(1.019),
csl::ComplexProperty::Complex
V_tb_mod->getValue(),
csl::ComplexProperty::Real
);
// As PDG says: angles in the first quadrant so cos > 0 and sin > 0
......@@ -357,7 +357,7 @@ namespace mty::sm_input {
inline
csl::Expr sin_CKM_13 = csl::Evaluated(
csl::sqrt_s(V_cb*V_cb + V_tb*V_tb),
csl::sqrt_s(1 - csl::pow_s(cos_CKM_13, 2)),
csl::eval::literal | csl::eval::numerical
);
......
......@@ -65,8 +65,14 @@ class FeynmanRule {
std::shared_ptr<wick::Graph>& getDiagram();
std::shared_ptr<wick::Graph> const& getDiagram() const;
csl::Expr getExpr() const;
bool contains(mty::QuantumFieldParent *parent) const;
size_t count(mty::QuantumFieldParent *parent) const;
bool contains(mty::QuantumFieldParent const *parent) const;
size_t count(mty::QuantumFieldParent const *parent) const;
bool contains(mty::Particle const &p) const {
return contains(p.get());
}
size_t count(mty::Particle const &p) const {
return count(p.get());
}
void setFieldProduct(std::vector<QuantumField> const& fieldProduct);
void setInteractionTerm(TermType const& term);
......@@ -103,6 +109,8 @@ class FeynmanRule {
csl::Expr expr;
};
int fermionicFactor(std::vector<mty::QuantumField> const &fieldProduct);
}
......
......@@ -64,6 +64,7 @@
#include "2HDM.h"
#include "MSSM.h"
#include "PMSSM.h"
#include "NMFV.h"
#include "wilson.h"
#include "wilsonOperator.h"
#include "propagator.h"
......
......@@ -221,6 +221,12 @@ protected:
void initInteractions172();
void initInteractions173();
void initInteractions174();
void initInteractions175();
void initInteractions176();
void initInteractions177();
void initInteractions178();
void initInteractions179();
void initInteractions180();
void initSpectrum();
protected:
......
......@@ -326,7 +326,7 @@ void PMSSM_LEM::initContent()
snu_tau->setMass(m_snu_tau);
addParticle(snu_tau, false);
st_1 = mty::scalarboson_s("st_1 ; \\tilde{st_1}", *this);
st_1 = mty::scalarboson_s("st_1 ; \\tilde{t_1}", *this);
st_1->setMass(m_st_L);
st_1->setGroupRep("C", {1, 0});
addParticle(st_1, false);
......@@ -336,7 +336,7 @@ void PMSSM_LEM::initContent()
st_2->setGroupRep("C", {1, 0});
addParticle(st_2, false);
sb_1 = mty::scalarboson_s("sb_1 ; \\tilde{sb_1}", *this);
sb_1 = mty::scalarboson_s("sb_1 ; \\tilde{b_1}", *this);
sb_1->setMass(m_sb_L);
sb_1->setGroupRep("C", {1, 0});
addParticle(sb_1, false);
......@@ -346,11 +346,11 @@ void PMSSM_LEM::initContent()
sb_2->setGroupRep("C", {1, 0});
addParticle(sb_2, false);
stau_1 = mty::scalarboson_s("stau_1 ; \\tilde{stau_1}", *this);
stau_1 = mty::scalarboson_s("stau_1 ; \\tilde{\\tau_1}", *this);
stau_1->setMass(m_stau_L);
addParticle(stau_1, false);
stau_2 = mty::scalarboson_s("stau_2 ; \\tilde{tau_2}", *this);
stau_2 = mty::scalarboson_s("stau_2 ; \\tilde{\\tau_2}", *this);
stau_2->setMass(m_stau_R);
addParticle(stau_2, false);
......@@ -931,6 +931,12 @@ void PMSSM_LEM::initInteractions()
initInteractions172();
initInteractions173();
initInteractions174();
initInteractions175();
initInteractions176();
initInteractions177();
initInteractions178();
initInteractions179();
initInteractions180();
}
void PMSSM_LEM::initSpectrum()
......
This diff is collapsed.
......@@ -14,9 +14,97 @@
// along with MARTY. If not, see <https://www.gnu.org/licenses/>.
#include <marty>
#include "pmssm_lem.h"
auto cc (csl::Expr e) { return csl::GetComplexConjugate(e); }
auto hc (csl::Expr e) { return csl::GetHermitianConjugate(e, &mty::dirac4); }
void print(csl::Expr e)
{
std::cout << "EXPR : " << e << '\n';
std::cout << "cc : " << cc(e) << '\n';
std::cout << "hc : " << hc(e) << '\n';
}
using namespace mty;
void square(Model & model, std::vector<Insertion> const &ins)
{
mty::option::displayAbbreviations = false;
auto ampl = model.computeAmplitude(TreeLevel, ins);
//for (auto &diag : ampl.getDiagrams())
// csl::Evaluate(diag.getExpression(), csl::eval::abbreviation);
std::cout << "FINAL AMPLITUDE : " << '\n';
Display(ampl);
//Show(ampl);
// auto square = model.computeSquaredAmplitude(ampl);
// std::cout << csl::Evaluated(square, csl::eval::abbreviation) << '\n';
DisplayAbbreviations("EXT");
}
int main() {
// PMSSM_LEM pmssm;
// std::cout << pmssm << '\n';
// auto amplmssm = pmssm.computeAmplitude(TreeLevel, {Incoming("C_1"), Incoming(AntiPart("N_1")), Outgoing("W")});
// for (auto &diag : amplmssm.getDiagrams())
// csl::Evaluate(diag.getExpression(), csl::eval::abbreviation);
// Display(amplmssm);
// Show(amplmssm);
// return 0;
using namespace csl;
Model model;
model.init();
auto N1 = diracfermion_s("N1", model);
auto N2 = diracfermion_s("N2", model);
auto N3 = diracfermion_s("N3", model);
N1->setSelfConjugate(true);
N2->setSelfConjugate(true);
auto C = diracfermion_s("C", model);
auto CL = C->getWeylFermion(Chirality::Left);
auto CR = C->getWeylFermion(Chirality::Right);
auto W = vectorboson_s("W", model);
//W->setSelfConjugate(false);
model.addParticles({N1, N2, N3, C, W});
auto gamma = dirac4.gamma;
auto Cm = dirac4.C_matrix;
auto a = DiracIndices(10);
auto mu = MinkowskiIndex();
auto aL = constant_s("a_L");
auto aR = constant_s("a_R");
model.addLagrangianTerm(
aL * W(mu) * cc(CR(a[0])) * gamma({+mu, a[0], a[1]}) * N2(a[1]),
true
);
model.addLagrangianTerm(
aR * W(mu) * gamma({+mu, a[1], a[2]}) * cc(N2(a[1]))*Cm({a[0],a[2]}) * cc(CL(a[0])),
true
);
//model.addLagrangianTerm(
// -aR * W(mu) * gamma({+mu, a[1], a[2]}) * cc(N2(a[1]))*Cm({a[2],a[0]}) * cc(CL(a[0])),
// true
// );
//model.addLagrangianTerm(
// aR * W(mu) *CR(a[0])* N2(a[2])*gamma({+mu,a[1],a[2]})*Cm({a[1],a[0]}),
// true
// );
model.refresh();
std::cout << model << '\n';
Display(model.getFeynmanRules());
auto ampl2 = model.computeAmplitude(TreeLevel, {Incoming("C"), Incoming("N2"), Outgoing("W")});
Display(ampl2);
//auto squared = model.computeSquaredAmplitude(ampl);
//std::cout << squared << '\n';
//std::cout << csl::Evaluated(squared, csl::eval::abbreviation) << '\n';
return 0;
///////////////////////////////////////////////////
// Setting the type of model we want, between
// 1 and 4 (inclusive)
......
......@@ -189,7 +189,7 @@ namespace sgl {
bool isTrace() const { return a == b && !(psiL || psiR); }
void checkGammaAndConjugation();
std::optional<GExpr> checkGammaAndConjugation() const;
static csl::Index easyIndex(int i) {
if (auto pos = m_easyIndex.find(i); pos != m_easyIndex.end())
......@@ -270,23 +270,22 @@ namespace sgl {
if (chain->isZero())
return cslexpr_s(CSL_0);
auto factor = chain->getFactor();
return (factor != CSL_1) ?
factor*chain->getTerm() :
GExpr(chain);
if (factor != CSL_1) {
auto term = std::dynamic_pointer_cast<IndexChain>(chain->getTerm());
auto opt_chain = term->checkGammaAndConjugation();
return factor * opt_chain.value_or(term);
}
return chain->checkGammaAndConjugation().value_or(chain);
}
template<class ...Args>
GExpr indexchain_s(
std::initializer_list<GExpr> gammas,
Args &&...args) {
auto chain = std::make_shared<IndexChain>(
std::vector<GExpr>{gammas}, std::forward<Args>(args)...);
if (chain->isZero())
return cslexpr_s(CSL_0);
auto factor = chain->getFactor();
return (factor != CSL_1) ?
factor*chain->getTerm() :
GExpr(chain);
return indexchain_s(
std::vector<GExpr>(gammas.begin(), gammas.end()),
std::forward<Args>(args)...
);
}
GExpr indexchain_s(
......
......@@ -54,7 +54,6 @@ namespace sgl {
a(t_a),
b(t_b)
{
checkGammaAndConjugation();
}
IndexChain::IndexChain(
......@@ -68,7 +67,6 @@ namespace sgl {
m_argument.reserve(mu.size());
for (const auto &i : mu)
m_argument.push_back(i);
checkGammaAndConjugation();
}
IndexChain::IndexChain(
......@@ -91,7 +89,6 @@ namespace sgl {
b(t_b),
psiL(field_s(t_a))
{
checkGammaAndConjugation();
}
IndexChain::IndexChain(
......@@ -106,7 +103,6 @@ namespace sgl {
m_argument.reserve(mu.size());
for (const auto &i : mu)
m_argument.push_back(i);
checkGammaAndConjugation();
}
IndexChain::IndexChain(
......@@ -129,7 +125,6 @@ namespace sgl {
b(t_b.index),
psiR(field_s(t_b))
{
checkGammaAndConjugation();
}
IndexChain::IndexChain(
......@@ -144,7 +139,6 @@ namespace sgl {
m_argument.reserve(mu.size());
for (const auto &i : mu)
m_argument.push_back(i);
checkGammaAndConjugation();
}
......@@ -170,7 +164,6 @@ namespace sgl {
psiL(field_s(t_a)),
psiR(field_s(t_b))
{
checkGammaAndConjugation();
}
IndexChain::IndexChain(
......@@ -186,11 +179,10 @@ namespace sgl {
m_argument.reserve(mu.size());
for (const auto &i : mu)
m_argument.push_back(i);
checkGammaAndConjugation();
}
void IndexChain::checkGammaAndConjugation()
std::optional<GExpr> IndexChain::checkGammaAndConjugation() const
{
for (const auto &arg : m_argument)
if (!IsType<GammaIndex>(arg)
......@@ -200,15 +192,22 @@ namespace sgl {
throw Exception::MathError;
}
if (m_argument.empty())
return;
return std::nullopt;
if (psiL && ConvertTo<GammaIndex>(m_argument.front())->isC()) {
psiL->conjugate();
m_argument.erase(m_argument.begin());
csl::Expr factor = psiL->isComplexConjugated() ? CSL_1 : CSL_M_1;
IndexChain cpy(*this);
cpy.psiL->conjugate();
cpy.m_argument.erase(cpy.m_argument.begin());
return cslexpr_s(factor)*cpy.copy();
}
if (psiR && ConvertTo<GammaIndex>(m_argument.back())->isC()) {
psiR->conjugate();
m_argument.erase(m_argument.end() - 1);
csl::Expr factor = !psiR->isComplexConjugated() ? CSL_1 : CSL_M_1;
IndexChain cpy(*this);
cpy.psiR->conjugate();
cpy.m_argument.erase(cpy.m_argument.end() - 1);
return cslexpr_s(factor)*cpy.copy();
}
return std::nullopt;
}
bool IndexChain::isZero() const
......@@ -1163,6 +1162,7 @@ namespace sgl {
sign *= -1;
}
});
LOG("RES :", sign, newChain.copy())
return { sign, newChain };
}
......
......@@ -1116,7 +1116,8 @@ void MSSM_Model::diagonalize2By2Matrices()
{eta_u, eta_d},
{G0, A0},
{{csl::sin_s(beta_h), csl::cos_s(beta_h)},
{-csl::cos_s(beta_h), csl::sin_s(beta_h)}}
{-csl::cos_s(beta_h), csl::sin_s(beta_h)}},
true, 1 // diagonalize, one massless state
);
mty::Particle phi_u = getParticle("phi_u");
......@@ -1127,7 +1128,8 @@ void MSSM_Model::diagonalize2By2Matrices()
{phi_u, phi_d},
{Gp, Hp},
{{csl::sin_s(beta_h), csl::cos_s(beta_h)},
{-csl::cos_s(beta_h), csl::sin_s(beta_h)}}
{-csl::cos_s(beta_h), csl::sin_s(beta_h)}},
true, 1 // diagonalize, one massless state