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

Color weights simplifications added

parent 9c742502
......@@ -551,6 +551,8 @@ class IndexStructure{
*/
std::vector<Index>::const_iterator find(const Index& t_index) const;
std::vector<Index>::iterator find(const Index& t_index);
/*!
* \return The part of the vector \b index that contains free indices.
*/
......
......@@ -445,6 +445,11 @@ vector<Index>::const_iterator IndexStructure::find(const Index& t_index) const
return std::find(begin(), end(), t_index);
}
vector<Index>::iterator IndexStructure::find(const Index& t_index)
{
return std::find(begin(), end(), t_index);
}
IndexStructure IndexStructure::getFreeStructure() const
{
// We create an empty structure and append only free indices to it
......
......@@ -1558,8 +1558,25 @@ bool TensorParent::hasContractionProperty(
const Abstract* self,
Expr_info B) const
{
if (IsIndicialTensor(B)
and B->getParent_info() == this
bool isBTensor = IsIndicialTensor(B);
csl::Parent_info BParent;
if (isBTensor)
BParent = B->getParent_info();
if (isBTensor) {
if ((BParent->getFullySymmetric() && getFullyAntiSymmetric())
|| (BParent->getFullyAntiSymmetric() && getFullySymmetric()))
{
size_t match = 0;
for (const auto &i : self->getIndexStructureView())
for (const auto &j : B->getIndexStructureView())
if (i == j)
++match;
if (match > 1)
return true;
}
}
if (isBTensor
and BParent == this
and not selfContraction.empty()) {
// We determine which contraction takes place between self and B
// and see if it appears in the list of contractions that have special
......@@ -1571,8 +1588,8 @@ bool TensorParent::hasContractionProperty(
return true;
}
}
else if (IsIndicialTensor(B)
and B->getParent_info() != this
else if (isBTensor
and BParent != this
and not extContraction.empty()) {
SelfContraction c(self, B);
for (size_t i = 0; i != extContraction.size(); ++i)
......@@ -1586,6 +1603,23 @@ bool TensorParent::hasContractionProperty(
Expr TensorParent::contraction(const Abstract* self, Expr_info B) const
{
bool isBTensor = IsIndicialTensor(B);
csl::Parent_info BParent;
if (isBTensor)
BParent = B->getParent_info();
if (isBTensor) {
if ((BParent->getFullySymmetric() && getFullyAntiSymmetric())
|| (BParent->getFullyAntiSymmetric() && getFullySymmetric()))
{
size_t match = 0;
for (const auto &i : self->getIndexStructureView())
for (const auto &j : B->getIndexStructureView())
if (i == j)
++match;
if (match > 1)
return CSL_0;
}
}
// This function should be called only if the contraction property has
// been verified with TensorParent::hasContractionProperty()
......
......@@ -494,6 +494,10 @@ namespace csl {
source << indent(2) << ");\n";
if (diag.squaredMass) {
for (const auto &mass : diag.masses) {
source << indent(1) << "if (0 > outputs." << mass << ") {\n";
source << indent(2) << "std::cerr << \"Warning: negative squared mass for \" << ";
source << '\"' << mass << '\"' << " << \".\\n\";\n";
source << indent(1) << "}\n";
source << indent(1) << "outputs." << mass << " = ";
source << ((quadruple) ? "sqrtq(" : "std::sqrt(");
source << "outputs." << mass << ");\n";
......
#pragma once
#include "colorSpace.h"
namespace mty {
csl::Expr CalculateFColorTrace(csl::Expr const &expr);
csl::Expr CalculateTrace(
const mty::SemiSimpleGroup* colorGroup,
const csl::Expr& expr);
csl::Expr CalculateColorTrace(csl::Expr const &init);
} // namespace mty
......@@ -184,13 +184,6 @@ class ColorSpace: public csl::Space {
const mty::SemiSimpleAlgebra* algebra;
};
csl::Expr CalculateTrace(const mty::SemiSimpleGroup* colorGroup,
const csl::Expr& expr);
} // End of namespace color
namespace mty {
csl::Expr CalculateColorTrace(csl::Expr const &init);
}
#endif
......@@ -220,6 +220,14 @@ class GaugedGroup {
mty::gauge::Choice getGaugeChoice() const;
csl::Expr getEffectiveQuadraticCasimir(
mty::Irrep const &irrep
) const;
csl::Expr getCA() const;
csl::Expr getCF() const;
mty::Generator getF() const;
csl::Tensor getD(size_t p) const;
......
......@@ -46,6 +46,7 @@
#include "diracology.h"
#include "CKM.h"
#include "colorSpace.h"
#include "colorSimplification.h"
#include "insertion.h"
#include "expander.h"
#include "gauge.h"
......
......@@ -17,6 +17,7 @@
#include "polarization.h"
#include "mrtInterface.h"
#include "colorSpace.h"
#include "colorSimplification.h"
#include "diracology.h"
#include "fermionChain.h"
#include "feynmanIntegral.h"
......
#include "colorSimplification.h"
#include "gaugedGroup.h"
namespace mty {
static std::set<color::ColorSpace const*> getColorSpaces(
csl::Expr const &expr
)
{
std::set<color::ColorSpace const*> colors;
csl::VisitEachLeaf(expr, [&](csl::Expr const &sub)
{
if (csl::IsIndicialTensor(sub)) {
std::vector<csl::Space const*> space = sub->getParent()->getSpace();
if (space.size() != 3
or space[0] == space[1]
or space[1] != space[2])
return;
auto const &pointed = *space[0];
if (typeid(color::ColorSpace) == typeid(pointed)
and typeid(color::ColorSpace) == typeid(pointed))
colors.insert(
static_cast<color::ColorSpace const*>(space[1]));
}
});
return colors;
}
static std::set<color::ColorSpace const*> getAdjointSpaces(
csl::Expr const &expr
)
{
std::set<color::ColorSpace const*> colors;
csl::VisitEachLeaf(expr, [&](csl::Expr const &sub)
{
if (csl::IsIndicialTensor(sub)) {
std::vector<csl::Space const*> space = sub->getParent()->getSpace();
if (space.size() != 3
or space[0] != space[1]
or space[1] != space[2])
return;
auto const &pointed = *space[0];
if (typeid(color::ColorSpace) == typeid(pointed)
and typeid(color::ColorSpace) == typeid(pointed))
colors.insert(
static_cast<color::ColorSpace const*>(space[1]));
}
});
return colors;
}
static std::vector<size_t> getPos_fABC(
csl::Expr &prod,
color::ColorSpace const *space
)
{
std::vector<size_t> pos_fABC;
const size_t sz = prod->size();
csl::Tensor f = space->getF();
pos_fABC.reserve(sz);
for (size_t i = 0; i != sz; ++i) {
if (csl::IsIndicialTensor(prod[i])
&& f.get() == prod[i]->getParent_info()) {
pos_fABC.push_back(i);
}
}
return pos_fABC;
}
static bool hasCommonIndex(
csl::Expr const &tensorA,
csl::Expr const &tensorB
)
{
for (const auto &i : tensorA->getIndexStructureView())
for (const auto &j : tensorB->getIndexStructureView())
if (i == j)
return true;
return false;
}
static bool isFTrace(
csl::Expr const &tensorA,
csl::Expr const &tensorB,
csl::Expr const &tensorC
)
{
return hasCommonIndex(tensorA, tensorB)
&& hasCommonIndex(tensorB, tensorC)
&& hasCommonIndex(tensorC, tensorA);
}
static int placeCommonIndex(
csl::IndexStructure &si,
csl::IndexStructure &sj,
size_t posIni,
size_t posInj
)
{
int sign = 1;
for (size_t i = 0; i != si.size(); ++i) { // index i1
if (auto posj = sj.find(si[i]); posj != sj.end()) {
if (i != posIni) {
sign *= -1;
std::swap(si[i], si[1]);
}
if (posj - sj.begin() != static_cast<long>(posInj)) {
sign *= -1;
std::swap(sj[0], *posj);
}
}
}
return sign;
}
static csl::Expr fTrace(
csl::Expr const &fi,
csl::Expr const &fj,
csl::Expr const &fk,
mty::GaugedGroup const *group
)
{
// Equation (44) of hep-ph/9802376v1
// Sorting indices as
// f_{i1 i2 a} * f_{i2 i3 b} * f_{i3 i1 c} to determine the sign
csl::IndexStructure si = fi->getIndexStructure();
csl::IndexStructure sj = fj->getIndexStructure();
csl::IndexStructure sk = fk->getIndexStructure();
int sign = 1;
sign *= placeCommonIndex(si, sj, 1, 0);
sign *= placeCommonIndex(sj, sk, 1, 0);
sign *= placeCommonIndex(sk, si, 1, 0);
csl::Index a = si[2];
csl::Index b = sj[2];
csl::Index c = sk[2];
csl::Parent f = fi->getParent();
csl::Expr CA = group->getCA();
return sign * CSL_HALF * CA * f({a, b, c});
}
static bool CalculateFColorTrace_impl(
csl::Expr &prod,
color::ColorSpace const *space
)
{
HEPAssert(csl::IsProd(prod),
mty::error::TypeError,
"Expected a product, " + toString(prod) + " given.")
std::vector<size_t> pos_fABC = getPos_fABC(prod, space);
if (pos_fABC.size() < 3)
return false;
for (size_t i = 0; i != pos_fABC.size(); ++i) {
auto &fi = prod[pos_fABC[i]];
for (size_t j = i+1; j < pos_fABC.size(); ++j) {
auto &fj = prod[pos_fABC[j]];
for (size_t k = j+1; k < pos_fABC.size(); ++k) {
auto &fk = prod[pos_fABC[k]];
if (isFTrace(fi, fj, fk)) {
fi = fTrace(fi, fj, fk, space->getGroup()->getGaugedGroup());
fj = fk = CSL_1;
return true;
}
}
}
}
return false;
}
csl::Expr CalculateFColorTrace(
csl::Expr const &expr
)
{
csl::Expr res = csl::DeepCopy(expr);
auto colors = getAdjointSpaces(expr);
csl::Transform(res, [&](csl::Expr &sub) {
if (csl::IsProd(sub)) {
for (const auto &c : colors)
return CalculateFColorTrace_impl(sub, c);
}
return false;
});
return res;
}
csl::Expr CalculateTrace(const mty::SemiSimpleGroup* colorGroup,
const csl::Expr& expr)
{
csl::ScopedProperty p(&csl::option::applyChainContractions, true);
std::vector<const csl::Space*> spaces = colorGroup->getAllVectorSpace();
for (const auto& space : spaces)
space->keepCycles = true;
csl::Expr res = DeepRefreshed(expr);
for (const auto& space : spaces)
space->keepCycles = false;
return res;
}
csl::Expr CalculateColorTrace(csl::Expr const &init)
{
csl::ScopedProperty p1(&csl::option::applyChainContractions, true);
auto colors = getColorSpaces(init);
for (auto &c : colors)
c->keepCycles = true;
csl::Expr res = csl::DeepRefreshed(init);
for (auto &c : colors)
c->keepCycles = false;
res = CalculateFColorTrace(res);
return res;
}
} // namespace mty
......@@ -719,50 +719,4 @@ tuple<csl::Index, csl::Index, int> ColorSpace::getFreeIndicesAdjoint(
return make_tuple(structA[0], structB[0], sign);
}
csl::Expr CalculateTrace(const mty::SemiSimpleGroup* colorGroup,
const csl::Expr& expr)
{
csl::ScopedProperty p(&csl::option::applyChainContractions, true);
vector<const csl::Space*> spaces = colorGroup->getAllVectorSpace();
for (const auto& space : spaces)
space->keepCycles = true;
csl::Expr res = DeepRefreshed(expr);
for (const auto& space : spaces)
space->keepCycles = false;
return res;
}
} // End of namespace color
namespace mty {
csl::Expr CalculateColorTrace(csl::Expr const &init)
{
csl::ScopedProperty p1(&csl::option::applyChainContractions, true);
std::set<color::ColorSpace const*> colors;
csl::VisitEachLeaf(init, [&](csl::Expr const &sub)
{
if (csl::IsIndicialTensor(sub)) {
std::vector<csl::Space const*> space = sub->getParent()->getSpace();
if (space.size() != 3
or space[0] == space[1]
or space[1] != space[2])
return;
auto const &pointed = *space[0];
if (typeid(color::ColorSpace) == typeid(pointed)
and typeid(color::ColorSpace) == typeid(pointed))
colors.insert(
static_cast<color::ColorSpace const*>(space[1]));
}
});
for (auto &c : colors)
c->keepCycles = true;
csl::Expr res = csl::DeepRefreshed(init);
for (auto &c : colors)
c->keepCycles = false;
return res;
}
}
......@@ -143,7 +143,16 @@ csl::Tensor GaugedGroup::getD(size_t p) const
std::vector<csl::Space const*> spaces(p, adjointSpace);
csl::Tensor d_p("d", spaces);
d_p->setFullySymmetric();
if (p == 3)
d_p->addTraceLessNess(spaces[0]); // d_AAB = 0 (~ Tr(T_B))
d[p] = d_p;
std::vector<csl::Index> indices(spaces.size());
for (auto &i : indices)
i = spaces[0]->generateIndex();
d_p->addSelfContraction(
d_p(indices), d_p(indices),
csl::constant_s("d_" + getName() + "_" + std::to_string(10*p + p))
);
return d_p;
}
......@@ -246,6 +255,28 @@ AlgebraState GaugedGroup::getHighestWeight(const csl::Space* space) const
return group->getHighestWeight(space);
}
csl::Expr GaugedGroup::getEffectiveQuadraticCasimir(mty::Irrep const &rep) const
{
auto defining = group->getAlgebra()->getDefiningRep();
csl::Expr res = normalization(getType())
* group->getAlgebra()->getQuadraticCasimir(rep.getHighestWeight())
/ group->getAlgebra()->getQuadraticCasimir(defining.getHighestWeight())
* rep.getDim() / defining.getDim();
return res;
}
csl::Expr GaugedGroup::getCA() const
{
Irrep adjoint = group->getAdjointRep();
return getEffectiveQuadraticCasimir(adjoint);
}
csl::Expr GaugedGroup::getCF() const
{
auto defining = group->getAlgebra()->getDefiningRep();
return getEffectiveQuadraticCasimir(defining);
}
csl::Expr GaugedGroup::covariantDerivative(const csl::Expr& field,
const Index& mu)
{
......@@ -460,11 +491,7 @@ void NonAbelianGauged::setGeneratorProperties(
csl::Tensor delta_adj = adjointSpace->getDelta();
csl::Tensor delta = space->getDelta();
// Quadratic Casimir
auto defining = group->getAlgebra()->getDefiningRep();
csl::Expr res = normalization(getType())
* group->getAlgebra()->getQuadraticCasimir(rep.getHighestWeight())
/ group->getAlgebra()->getQuadraticCasimir(defining.getHighestWeight())
* rep.getDim() / defining.getDim();
csl::Expr res = getEffectiveQuadraticCasimir(rep);
if (space == adjointSpace)
res *= -1;
T->addSelfContraction(
......@@ -539,11 +566,8 @@ mty::gauge::GroupType SUGauged::getType() const
void SUGauged::init_f_d_ABC()
{
csl::Tensor d_ABC = csl::tensor_s("d_"+group->getName(),
{adjointSpace,
adjointSpace,
adjointSpace});
d_ABC->setFullySymmetric();
csl::Tensor d_ABC = getD(3);
d_ABC->addTraceLessNess(adjointSpace);
const size_t N = group->getAlgebra()->getOrderL() + 1;
if (N == 2) {
......@@ -569,7 +593,6 @@ void SUGauged::init_f_d_ABC()
f->setTensor(gell_mann::f);
d_ABC->setTensor(gell_mann::d3);
}
d[3] = d_ABC;
}
void SUGauged::setGeneratorProperties(
......
......@@ -14,6 +14,7 @@
// along with MARTY. If not, see <https://www.gnu.org/licenses/>.
#include "quantumField.h"
#include "ghostField.h"
#include "kinematics.h"
#include "insertion.h"
#include "mrtError.h"
......@@ -81,7 +82,8 @@ namespace mty {
return A.getField() == B.getField()
&& A.isIncoming() == B.isIncoming()
&& (A.isParticle() == B.isParticle()
|| A.getField()->isSelfConjugate());
|| A.getField()->isSelfConjugate()
|| IsOfType<GhostBoson>(A.getField()));
}
csl::Expr Kinematics::getDegreesOfFreedomFactor() const
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment