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

Delta F = 2 op, sign and factor Ok, sqrt(2) issue fixed

parent de5345e8
......@@ -1794,6 +1794,10 @@ bool Prod::mergeTerms()
matched = false;
continue;
}
if (csl::IsNumerical(argument[i]) || csl::IsNumerical(argument[i+1])) {
matched = false;
continue;
}
if (not matched) {
factor = CSL_1;
if (argument[i]->getType() == csl::Type::Pow) { //Pow
......@@ -1828,9 +1832,9 @@ bool Prod::mergeTerms()
if (*term==term2) {
matched = true;
factor = factor->addition_own(factor2); //factor->addition_own(factor2);
if (term->isInteger()
and factor->getType() == csl::Type::IntFraction)
continue;
//if (term->isInteger()
// and factor->getType() == csl::Type::IntFraction)
// continue;
argument.erase(argument.begin()+i+1);
argument[i] = pow_s(term, factor);
if (*argument[i] == CSL_1) {
......
......@@ -310,7 +310,9 @@ namespace mty {
*
* @param t_kinematics New kinematic context for the process.
*/
void setKinematics(Kinematics const &t_kinematics);
void setKinematics(
Kinematics const &t_kinematics,
bool replace = true);
private:
......
......@@ -125,9 +125,10 @@ namespace mty {
* @param kinematics Kinematics that will follow the same order as the
* insertions.
*/
void orderInsertions(
static void orderInsertions(
std::vector<mty::QuantumField> &insertions,
Kinematics &kinematics
Kinematics &kinematics,
FeynOptions &options
);
static void applyMomentumVertices(
......
......@@ -176,6 +176,11 @@ class FeynmanIntegral: public csl::AbstractMultiFunc {
std::vector<csl::Expr> const& masses,
std::vector<csl::Index> const& indices);
static void removeExternalMomenta(
csl::Expr &expr,
csl::Parent const &Q
);
FeynmanIntegral(IntegralType t_type,
int t_looptoolsId,
std::vector<csl::Expr> const& t_argument);
......
......@@ -50,6 +50,12 @@ class FeynOptions;
class Kinematics;
class Amplitude;
enum class DecompositionMode {
Minimal,
BasisProjection,
Matching
};
/*!
* \brief Contains all objects in the theory.
* In particular QuantumField objects, Gauge, Flavor, Particle...
......@@ -166,7 +172,7 @@ public:
std::vector<Lagrangian::TermType> &lagrangian,
std::vector<Insertion> insertions,
Kinematics const &kinematics,
FeynOptions const &options,
FeynOptions options,
std::vector<FeynmanRule const*> rules = {}
);
......@@ -202,18 +208,18 @@ public:
WilsonSet getWilsonCoefficients(
Amplitude const &ampl,
FeynOptions const &feynOptions,
bool squaredAfter = false
DecompositionMode mode = DecompositionMode::Matching
);
WilsonSet getWilsonCoefficients(
Amplitude const &ampl,
bool squaredAfter = false
DecompositionMode mode = DecompositionMode::Matching
);
WilsonSet computeWilsonCoefficients(
int order,
std::vector<Insertion> const &insertions,
FeynOptions const &feynOptions = {}
FeynOptions feynOptions = {}
);
WilsonSet computeWilsonCoefficients_default(
......@@ -246,8 +252,8 @@ public:
);
WilsonSet computeWilsonCoefficients_4Fermions(
std::vector<Insertion> const &insertions,
FeynOptions const &feynOptions = {}
std::vector<Insertion> insertions,
FeynOptions feynOptions = {}
);
Amplitude connectAmplitudes(
......@@ -321,6 +327,13 @@ protected:
int operatorDegeneracy(std::vector<mty::Insertion> const &insertions);
int matchingFermionSign(std::vector<int> fermionOrder);
int fermionSign(
std::vector<Insertion> const &model,
std::vector<Insertion> order
);
} // End of namespace mty
#endif /* MODE_H_INCLUDED */
......@@ -529,6 +529,7 @@ void Display(WilsonSet const& wilsons,
void Show(std::vector<FeynmanRule> const& rules);
void Show(std::vector<std::shared_ptr<wick::Graph>> const& diagrams);
void Show(mty::Amplitude const& diagrams);
void Show(WilsonSet const &wilsons);
void Show(
std::vector<FeynmanRule> const& rules,
......@@ -545,6 +546,11 @@ void Show(
size_t first,
size_t last
);
void Show(
WilsonSet const &wilsons,
size_t first,
size_t last
);
template<class T>
void Display(T const& printable,
......
......@@ -142,6 +142,7 @@ struct WilsonSet: public std::vector<Wilson> {
Kinematics kinematics {};
FeynOptions options {};
std::vector<std::shared_ptr<wick::Graph>> graphs;
};
std::vector<Wilson> copyWilsons(std::vector<Wilson> const &wilsons);
......
......@@ -168,12 +168,13 @@ namespace sgl {
csl::Index a = expr->getIndexStructureView()[nI-2];
csl::Index b = expr->getIndexStructureView()[nI-1];
minkoIndices.erase(minkoIndices.end() - 2, minkoIndices.end());
LOG(indexchain_s(
csl::Expr factor = (minkoIndices.size() == 2) ? CSL_I : CSL_1;
LOG(factor*indexchain_s(
gammaindex_s(minkoIndices),
a,
b
))
return indexchain_s(
return factor*indexchain_s(
gammaindex_s(minkoIndices),
a,
b
......
......@@ -152,9 +152,10 @@ namespace sgl {
}
csl::Tensor gamma = pos->second;
std::vector<csl::Index> indices(m_indices);
csl::Expr factor = (indices.size() == 2) ? -CSL_I : CSL_1;
indices.push_back(a);
indices.push_back(b);
return gamma(indices);
return factor*gamma(indices);
}
csl::Expr GammaIndex::toCSL(TensorSet const &) const
......
......@@ -860,7 +860,7 @@ namespace sgl {
};
case 4:
return {
cslexpr_s(CSL_HALF),
cslexpr_s(-CSL_HALF),
IndexChain(gammaindex_s({mu, nu}), e[0], e[1]),
IndexChain(gammaindex_s({+mu, +nu}), e[2], e[3])
};
......@@ -903,7 +903,7 @@ namespace sgl {
};
case 4:
return {
cslexpr_s(1),
cslexpr_s(-1),
IndexChain(gammaindex_s({mu, nu}), e[0], e[1]),
IndexChain(gammaindex_s({+mu, +nu}), e[2], e[3])
};
......
......@@ -153,10 +153,12 @@ namespace mty {
return Amplitude { options, diagramCopy, kinematics };
}
void Amplitude::setKinematics(Kinematics const &t_kinematics)
void Amplitude::setKinematics(Kinematics const &t_kinematics, bool replace)
{
for (auto &diag : diagrams) {
Kinematics::replace(diag.getExpression(), kinematics, t_kinematics);
if (replace) {
for (auto &diag : diagrams) {
Kinematics::replace(diag.getExpression(), kinematics, t_kinematics);
}
}
kinematics = t_kinematics;
}
......
......@@ -94,7 +94,7 @@ namespace mty {
true
);
if (options.orderInsertions)
orderInsertions(insertions, kinematics);
orderInsertions(insertions, kinematics, options);
csl::PseudoIntegral integral;
integral.addTerm(LSZInsertion);
......@@ -160,7 +160,7 @@ namespace mty {
false);
// std::cout << "LSZ = " << LSZInsertion << '\n';
if (options.orderInsertions)
orderInsertions(insertions, kinematics);
orderInsertions(insertions, kinematics, options);
csl::PseudoIntegral integral;
integral.addTerm(LSZInsertion);
......@@ -429,7 +429,8 @@ namespace mty {
void AmplitudeInitializer::orderInsertions(
std::vector<mty::QuantumField> &insertions,
Kinematics &kinematics
Kinematics &kinematics,
FeynOptions &options
)
{
std::vector<size_t> indices(kinematics.size());
......
......@@ -704,7 +704,8 @@ namespace mty::simpli {
};
csl::DeepPartialExpand(expr, isEmitter, isReceiver);
csl::DeepFactor(expr);
// Do not add factorization here as the epsilon is not simplified
// right away and it could simply refactor what has been expanded
return (init != expr);
}
......@@ -1033,10 +1034,10 @@ namespace mty::simpli {
{
csl::ScopedProperty commut(&csl::option::checkCommutations, false);
csl::Refresh(expr);
if (mode != WilsonCoefficient && mode != FeynmanRule) {
//MomentumConservater cons(insertions, momenta);
//cons.apply(expr);
}
//if (mode != WilsonCoefficient && mode != FeynmanRule) {
// MomentumConservater cons(insertions, momenta);
// cons.apply(expr);
//}
std::vector<csl::Expr> factors;
if (csl::IsProd(expr)) {
for (size_t i = 0; i != expr->size(); ++i) {
......@@ -1064,8 +1065,6 @@ namespace mty::simpli {
if (mode == Amplitude)
reduceTensorIntegrals(expr);
//if (mode != SquaredAmplitude && mode != FeynmanRule)
// simplifyImpulsions(expr, insertions, momenta);
expandMomentaExperimental(expr, momenta);
size_t maxLoop = 10;
bool simplified;
......@@ -1110,8 +1109,6 @@ namespace mty::simpli {
addLocalTerms(expr);
csl::DeepFactor(expr);
}
// std::cout << "Before Abbrevs" << '\n';
// std::cout << expr << '\n';
//////////////////////
......@@ -1139,7 +1136,6 @@ namespace mty::simpli {
expr *= factor;
}
csl::DeepRefresh(expr);
// std::cout << "Factors = " << prod_s(factors) << '\n';
if (mode == Amplitude) {
abbreviateIntegral(expr);
......@@ -1147,8 +1143,6 @@ namespace mty::simpli {
abbreviateAll(expr);
csl::DeepRefresh(expr);
}
// std::cout << "RES = " << '\n';
// std::cout << expr << '\n';
}
} // namespace mty
......@@ -540,9 +540,11 @@ namespace mty {
std::string const &group
)
{
for (const auto &coupling : colorCouplings)
if (coupling.groupName == group)
for (const auto &coupling : colorCouplings) {
if (coupling.groupName == group) {
return coupling.coupling;
}
}
return ColorCoupling::Id;
}
......@@ -641,6 +643,8 @@ namespace mty {
};
if (!requirementsSatisfied(insertions, requirements))
return {};
for (const auto &color : colorCouplings)
model.getGroup(color.groupName); // Checking that the color exists
setUpDim6FermionOrder(wilsons, fermionOrder);
auto [psi1, psi2, psi3, psi4] = getExternalDim6Fermions(
insertions, fermionOrder);
......
......@@ -964,8 +964,8 @@ csl::Expr FeynmanIntegral::replaceIntegral(csl::Expr const& expr)
+ " given.");
csl::Parent variable = integral->getParent();
csl::Expr res = csl::DeepExpandedIf(integral->getOperand(),
[&](csl::Expr const& el)
csl::Expr res = integral->getOperand();
csl::DeepExpandIf(res, [&](csl::Expr const& el)
{
return el->dependsExplicitlyOn(variable.get());
});
......@@ -1049,6 +1049,14 @@ csl::Expr FeynmanIntegral::replaceIntegral(
&& static_cast<int>(indices.size()) < Ncrit) {
// Also check that non divergent integral
std::fill(momentum.begin(), momentum.end(), CSL_0);
// Removing external momenta in the numerator also
for (auto &f : factor) {
removeExternalMomenta(f, variable);
csl::DeepRefresh(f);
if (f == CSL_0) {
return CSL_0;
}
}
}
size_t firstZero = size_t(-1);
......@@ -1196,6 +1204,24 @@ csl::Expr FeynmanIntegral::applyQSquared(
return term1 + term2;
}
void FeynmanIntegral::removeExternalMomenta(
csl::Expr &expr,
csl::Parent const &Q
)
{
if (csl::IsSum(expr) || csl::IsProd(expr) || csl::IsPow(expr)) {
for (size_t i = 0; i != expr->size(); ++i)
removeExternalMomenta(expr[i], Q);
}
else if (csl::IsIndicialTensor(expr)) {
auto const &parent = expr->getParent();
auto const &spaces = parent->getSpace();
if (spaces.size() == 1 && spaces[0] == &csl::Minkowski && parent != Q) {
expr = CSL_0;
}
}
}
void testError(IntegralType type,
std::vector<csl::Expr> const& momentum,
size_t nMomentum,
......
......@@ -128,7 +128,8 @@ FeynmanRule::FeynmanRule(
fieldParent->setEnabledInDiagrams(false);
return;
}
auto wilsons = model.getWilsonCoefficients(res, options, true);
auto wilsons = model.getWilsonCoefficients(
res, options, DecompositionMode::Minimal);
diagram = res.getDiagrams()[0].getDiagram();
std::vector<csl::Expr> expressions = wilsons.obtainExpressions();
for (auto& ampl : expressions) {
......
......@@ -150,13 +150,14 @@ mty::Amplitude Model::computeAmplitude(
std::vector<Lagrangian::TermType> &lagrangian,
std::vector<mty::Insertion> insertions,
Kinematics const &kinematics,
FeynOptions const &options,
FeynOptions options,
std::vector<FeynmanRule const*> rules
)
{
std::vector<int> fermionOrder = options.getFermionOrder();
if (fermionOrder.empty())
fermionOrder = defaultFermionOrder(insertions);
options.setFermionOrder(fermionOrder);
if (options.orderExternalFermions)
applyFermionOrder(insertions, fermionOrder);
auto quantumInsertions = recoverQuantumInsertions(GetExpression(insertions));
......@@ -248,7 +249,8 @@ csl::Expr Model::computeSquaredAmplitude(
FeynOptions optionsL = amplL.getOptions();
// optionsL.setWilsonOperatorBasis(OperatorBasis::None);
auto wilsonsL = getWilsonCoefficients(amplL, optionsL, true);
auto wilsonsL = getWilsonCoefficients(
amplL, optionsL, DecompositionMode::Minimal);
if (&amplL == &amplR) { // same ampl, no need to recalculate
return computeSquaredAmplitude(
......@@ -259,7 +261,8 @@ csl::Expr Model::computeSquaredAmplitude(
}
FeynOptions optionsR = amplR.getOptions();
// optionsR.setWilsonOperatorBasis(OperatorBasis::None);
auto wilsonsR = getWilsonCoefficients(amplR, optionsR, true);
auto wilsonsR = getWilsonCoefficients(
amplR, optionsR, DecompositionMode::Minimal);
return computeSquaredAmplitude(
wilsonsL,
......@@ -397,10 +400,10 @@ void Model::projectOnBasis(
WilsonSet Model::getWilsonCoefficients(
Amplitude const &ampl,
bool squaredAfter
DecompositionMode mode
)
{
return getWilsonCoefficients(ampl, ampl.getOptions(), squaredAfter);
return getWilsonCoefficients(ampl, ampl.getOptions(), mode);
}
int operatorDegeneracy(std::vector<mty::Insertion> const &insertions)
......@@ -424,21 +427,42 @@ int operatorDegeneracy(std::vector<mty::Insertion> const &insertions)
return factor;
}
static int nPermutations(std::vector<int> &fermionOrder)
{
for (size_t i = 0; i+1 < fermionOrder.size(); ++i) {
if (fermionOrder[i+1] < fermionOrder[i]) {
std::swap(fermionOrder[i], fermionOrder[i+1]);
return 1 + nPermutations(fermionOrder);
}
}
return 0;
}
int matchingFermionSign(std::vector<int> fermionOrder)
{
return (nPermutations(fermionOrder) & 1) ? -1 : 1;
}
WilsonSet Model::getWilsonCoefficients(
Amplitude const &ampl,
FeynOptions const &feynOptions,
bool squaredAfter
DecompositionMode mode
)
{
mty::option::displayAbbreviations = false;
mty::cachedWilson.clear();
csl::ScopedProperty commut(&csl::option::checkCommutations, false);
std::vector<csl::Expr> amplitudesfull(ampl.obtainExpressions());
const auto basis = OperatorBasis::Standard;
const auto factor = feynOptions.getWilsonOperatorCoefficient();
csl::Abbrev::enableGenericEvaluation("EXT");
bool isMinimal = (mode == DecompositionMode::Minimal);
for (auto &ampl : amplitudesfull) {
if (!squaredAfter) {
ampl = CSL_I * csl::Copy(ampl);
ampl = csl::DeepCopy(ampl);
if (mode == DecompositionMode::Matching) {
ampl = CSL_I * ampl;
}
if (!isMinimal) {
projectOnBasis(ampl, basis);
}
csl::Evaluate(ampl);
......@@ -451,14 +475,12 @@ WilsonSet Model::getWilsonCoefficients(
auto wilsons = match(
amplitudesfull,
factor,
!squaredAfter && (basis == OperatorBasis::Standard),
squaredAfter);
!isMinimal && (basis == OperatorBasis::Standard),
isMinimal);
auto const &insertions = ampl.getKinematics().getInsertions();
auto const &momenta = ampl.getKinematics().getMomenta();
const int degeneracy = operatorDegeneracy(insertions);
if (!feynOptions.getFeynRuleCalculation()) {
std::cout << "Simplifying wilsons ..." << std::endl;
}
const int effSign = matchingFermionSign(feynOptions.getFermionOrder());
csl::ProgressBar bar(wilsons.size());
size_t index = 0;
std::vector<csl::Expr> coefs;
......@@ -468,7 +490,7 @@ WilsonSet Model::getWilsonCoefficients(
bar.progress(index++);
}
csl::Expr a = w.coef.getCoefficient();
if (!squaredAfter)
if (!isMinimal)
csl::Replace(a, csl::DMinko, csl::int_s(4));
if (amplitudesfull.size() < 2000) {
// Only if less than 2000 amplitudes
......@@ -477,7 +499,7 @@ WilsonSet Model::getWilsonCoefficients(
recoverQuantumInsertions(GetExpression(insertions)),
momenta,
ampl.getOptions(),
(squaredAfter) ?
(isMinimal) ?
simpli::SquaredAmplitude : simpli::WilsonCoefficient
);
csl::DeepHardFactor(a);
......@@ -489,13 +511,13 @@ WilsonSet Model::getWilsonCoefficients(
csl::Expr &a = coefs[i];
Wilson &w = wilsons[i];
csl::matcher::compress(a, 2);
if (squaredAfter) {
if (isMinimal) {
a = csl::Abbrev::makeAbbreviation("Cw", a);
}
else {
else if (mode == DecompositionMode::Matching) {
// Taking into account degeneracy if it is a real Wilson coefficient
// calculation
a /= csl::int_s(degeneracy);
a *= csl::intfraction_s(effSign, degeneracy);
}
w.coef.setCoefficient(a);
}
......@@ -504,6 +526,7 @@ WilsonSet Model::getWilsonCoefficients(
//res.mergeSorted();
wilsons.options = feynOptions;
wilsons.kinematics = ampl.getKinematics();
wilsons.graphs = ampl.obtainGraphs();
return wilsons;
}
......@@ -511,7 +534,7 @@ WilsonSet Model::getWilsonCoefficients(
WilsonSet Model::computeWilsonCoefficients(
int order,
std::vector<Insertion> const &insertions,
FeynOptions const &feynOptions
FeynOptions feynOptions
)
{
if (order == TreeLevel) {
......@@ -535,10 +558,16 @@ WilsonSet Model::computeWilsonCoefficients(
return computeWilsonCoefficients_2Fermions_1Vector(
insertions, feynOptions
);
if (nTot == 4 && nF == 4)
if (nTot == 4 && nF == 4) {
if (feynOptions.getFermionOrder().empty()) {
CallHEPError(mty::error::RuntimeError,
"A fermion order order must be provided for the dim 6 Wilson "
"coefficient calculation through the FeynOptions object.")
}
return computeWilsonCoefficients_4Fermions(
insertions, feynOptions
);
}
return computeWilsonCoefficients_default(
OneLoop, insertions, feynOptions
);
......@@ -574,12 +603,35 @@ WilsonSet Model::computeWilsonBoxes_4Fermions(
// feynOptions.verboseAmplitude = false;
feynOptions.setTopology(feynOptions.getTopology() & Topology::Box);
auto ampl = computeAmplitude(
OneLoop, kinematics.getInsertions(), feynOptions);
ampl.setKinematics(kinematics.alignedWith(ampl.getKinematics()));
auto wilsons = getWilsonCoefficients(ampl, feynOptions);
OneLoop, kinematics.getInsertions(), kinematics, feynOptions);
// ampl.setKinematics(kinematics.alignedWith(ampl.getKinematics()));
auto wilsons = getWilsonCoefficients(ampl); //, feynOptions);
return wilsons;
}
int fermionSign(
std::vector<Insertion> const &model,
std::vector<Insertion> order
)
{
// Supposes only fermion insertions with same size model && order
int nPerm = 0;
for (size_t i = 0; i != model.size(); ++i) {
size_t index = i;
for (size_t j = i+1; j != order.size(); ++j) {
if (model[i] == order[j]) {
index = j;
break;
}
}
if (index != i) {
std::swap(order[i], order[index]);
nPerm += index - i;
}
}
return (nPerm & 1) ? -1 : 1;