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

Wilson coefficient interface 1.0 !

parent 7e6fe4bf
......@@ -139,11 +139,11 @@ testDependencies_MACOS()
{
echo "Testing MACOS dependencies"
vmac=${v:1:${#v}-1}
testDependency_MACOS coreutils
testDependency_MACOS gcc@$vmac
testDependency_MACOS gsl
testDependency_MACOS qt5
testDependency_MACOS qt-creator cask
testDependency_MACOS coreutils
testDependency_MACOS gsl
message "${GREEN}${BOLD}[1 / 5] All dependencies are ready for compilation${NC}${NORMAL}"
}
......
......@@ -446,13 +446,7 @@ Expr LibEvalSession::addEval(Expr init)
}
auto pos = std::find(eval.begin(), eval.end(), term);
if (pos != eval.end()) {
//Expr factor = init->getNumericalFactor()
// / pos->init->getNumericalFactor();
//return factor * DeepCopy(pos->getExpr());
return Copy(pos->getExpr());
// init->getNumericalFactor()
// / pos->init->getNumericalFactor()
// * Copy(pos->getExpr());
}
csl::DeepRefresh(term);
LibEval newEval(term, maxID++);
......
......@@ -304,7 +304,7 @@ namespace csl {
printNameMapElement(
out, nIndent, param, LibraryGenerator::realUsing);
}
out << "};\n\n";
out << LibraryGenerator::indent(nIndent) << "};\n\n";
out << LibraryGenerator::indent(nIndent)
<< "std::map<std::string, csl::InitSanitizer<"
<< LibraryGenerator::complexUsing << ">*> complexParams {\n";
......
......@@ -125,7 +125,7 @@ namespace mty {
* @param kinematics Kinematics that will follow the same order as the
* insertions.
*/
static void orderInsertions(
void orderInsertions(
std::vector<mty::QuantumField> &insertions,
Kinematics &kinematics
);
......@@ -152,7 +152,7 @@ namespace mty {
inline static std::vector<mty::FeynmanRule const*> defaultFR{};
mty::Model const *model;
mty::FeynOptions const &options;
mty::FeynOptions options;
std::vector<mty::QuantumField> insertions;
Kinematics kinematics;
......
#pragma once
#include <string>
#include <vector>
namespace csl {
class Expr;
class Index;
}
namespace mty {
class Model;
struct Wilson;
struct WilsonSet;
enum class DiracCoupling {
S, P, L, R,
V, A, VL, VR,
T, TA, TL, TR
};
enum class ColorCoupling {
Id,
Generator,
Crossed,
InvCrossed // relevant only with conjugate representations (not SM)
};
struct ColorSpec {
std::string groupName;
ColorCoupling coupling;
};
std::vector<Wilson> getMagneticCoupling(
DiracCoupling coupling,
csl::Index alpha,
csl::Index beta
);
csl::Expr getMagneticGenerator(
mty::Model const &model,
csl::Expr &psi_star,
csl::Expr &psi,
csl::Expr &A
);
std::vector<Wilson> chromoMagneticOperator(
Model const &model,
WilsonSet const &wilsons,
DiracCoupling coupling
);
std::vector<Wilson> dimension6Operator(
Model const &model,
WilsonSet const &wilsons,
DiracCoupling leftCurrent,
DiracCoupling rightCurrent,
std::vector<ColorSpec> const &colorCouplings,
std::vector<int> fermionOrder = {}
);
inline std::vector<Wilson> dimension6Operator(
Model const &model,
WilsonSet const &wilsons,
DiracCoupling leftCurrent,
DiracCoupling rightCurrent,
std::vector<int> fermionOrder = {}
)
{
return dimension6Operator(
model,
wilsons,
leftCurrent,
rightCurrent,
std::vector<ColorSpec>{},
fermionOrder
);
}
inline std::vector<Wilson> dimension6Operator(
Model const &model,
WilsonSet const &wilsons,
DiracCoupling leftCurrent,
DiracCoupling rightCurrent,
ColorSpec const &colorCoupling,
std::vector<int> fermionOrder = {}
)
{
return dimension6Operator(
model,
wilsons,
leftCurrent,
rightCurrent,
std::vector<ColorSpec>{colorCoupling},
fermionOrder
);
}
}
......@@ -36,16 +36,20 @@ namespace mty {
mty::Insertion const &right
);
std::vector<mty::Insertion const*> fermionsOf(
std::vector<mty::Insertion> const &fields
);
std::vector<mty::Insertion*> fermionsOf(
std::vector<mty::Insertion> &fields
);
std::vector<int> defaultFermionOrder(
std::vector<mty::Insertion*> const &fields
std::vector<mty::Insertion const*> const &fields
);
inline std::vector<int> defaultFermionOrder(
std::vector<mty::Insertion> &fields
std::vector<mty::Insertion> const &fields
)
{
return defaultFermionOrder(fermionsOf(fields));
......
......@@ -67,6 +67,7 @@
#include "NMFV.h"
#include "wilson.h"
#include "wilsonOperator.h"
#include "builtinOperators.h"
#include "propagator.h"
#include "mrtUtils.h"
#include "wick.h"
......
......@@ -111,14 +111,25 @@ public:
mty::Amplitude computeAmplitude(
int order,
std::vector<mty::Insertion> insertions,
FeynOptions options = {}
FeynOptions &options
);
mty::Amplitude computeAmplitude(
int order,
std::vector<mty::Insertion> insertions,
Kinematics const &kinematics,
FeynOptions options = {}
FeynOptions &options
);
mty::Amplitude computeAmplitude(
int order,
std::vector<mty::Insertion> insertions
);
mty::Amplitude computeAmplitude(
int order,
std::vector<mty::Insertion> insertions,
Kinematics const &kinematics
);
mty::Amplitude computePartialAmplitude(
......
#pragma once
#include <utility>
namespace csl {
class Expr;
}
namespace mty {
std::pair<csl::Expr, csl::Expr> calculate_C9_C10();
}
......@@ -122,6 +122,11 @@ struct WilsonSet: public std::vector<Wilson> {
using std::vector<Wilson>::vector;
WilsonSet(WilsonSet const &) = default;
WilsonSet(WilsonSet &&) = default;
WilsonSet &operator=(WilsonSet const &) = default;
WilsonSet &operator=(WilsonSet &&) = default;
std::vector<csl::Expr> obtainExpressions() const {
std::vector<csl::Expr> res(size());
for (size_t i = 0; i != size(); ++i) {
......
......@@ -26,7 +26,24 @@
#include "wilson.h"
#include "diracology.h"
#include "polarization.h"
#include "builtinOperators.h"
namespace mty {
csl::Expr getWilsonCoefficient(
WilsonSet const &wilsons,
std::vector<Wilson> const &contributions
);
csl::Expr getWilsonCoefficient(
Model const &model,
WilsonSet const &wilsons,
DiracCoupling coupling
);
} // End of namespace mty
// Old way
namespace mty {
enum class OperatorType {
......@@ -136,6 +153,6 @@ namespace OperatorParser {
} // End of namespace OperatorParser
} // End of namespace mty
}
#endif /* ifndef MTY_WILSONOPERATOR_H_INCLUDED */
......@@ -38,8 +38,7 @@ namespace mty {
lagrangian(t_lagrangian),
feynmanRules(defaultFR)
{
}
}
AmplitudeInitializer::AmplitudeInitializer(
mty::Model const *t_model,
......@@ -435,6 +434,7 @@ namespace mty {
{
std::vector<size_t> indices(kinematics.size());
std::iota(begin(indices), end(indices), 0);
std::vector<int> fermionOrder = options.getFermionOrder();
for (size_t i = 0; i != insertions.size(); ++i) {
size_t mini = i;
for (size_t j = i+1; j < insertions.size(); ++j) {
......@@ -444,12 +444,16 @@ namespace mty {
if (mini != i) {
std::swap(insertions[i], insertions[mini]);
std::swap(indices[i], indices[mini]);
for (int &pos : fermionOrder) {
if (pos == static_cast<int>(mini)) pos = i;
else if (pos == static_cast<int>(i)) pos = mini;
}
}
}
kinematics.applyPermutation(indices);
options.setFermionOrder(fermionOrder);
}
void AmplitudeInitializer::initMomentumVertices(
std::vector<FeynmanRule> &localRules,
std::map<csl::Tensor, size_t> &vertexIds,
......
This diff is collapsed.
......@@ -43,6 +43,20 @@ namespace mty {
return false;
}
std::vector<mty::Insertion const*> fermionsOf(
std::vector<mty::Insertion> const &fields
)
{
std::vector<mty::Insertion const*> qfields;
qfields.reserve(fields.size());
for (auto &f : fields) {
if (f.getField()->isFermionic())
qfields.push_back(&f);
}
return qfields;
}
std::vector<mty::Insertion*> fermionsOf(
std::vector<mty::Insertion> &fields
)
......@@ -60,7 +74,7 @@ namespace mty {
void applyGoodAdjacentOrder(
std::vector<int> &numbers,
std::vector<size_t> &indicesLeft,
std::vector<mty::Insertion*> const &fermions
std::vector<mty::Insertion const*> const &fermions
)
{
for (size_t i = 0; i != indicesLeft.size(); i += 2) {
......@@ -83,7 +97,7 @@ namespace mty {
void applyGoodNonAdjacentOrder(
std::vector<int> &numbers,
std::vector<size_t> &indicesLeft,
std::vector<mty::Insertion*> const &fermions
std::vector<mty::Insertion const*> const &fermions
)
{
for (size_t i = 0; i != indicesLeft.size(); i += 2) {
......@@ -113,7 +127,7 @@ namespace mty {
}
std::vector<int> defaultFermionOrder(
std::vector<mty::Insertion*> const &fields
std::vector<mty::Insertion const*> const &fields
)
{
std::vector<int> numbers(fields.size());
......@@ -158,7 +172,7 @@ namespace mty {
{
applyFermionOrder(
insertions,
defaultFermionOrder(fermionsOf(insertions))
defaultFermionOrder(insertions)
);
}
}
......@@ -86,21 +86,42 @@ std::vector<mty::QuantumField> Model::recoverQuantumInsertions(
return fields;
}
mty::Amplitude Model::computeAmplitude(
int order,
std::vector<mty::Insertion> insertions
)
{
FeynOptions options;
return computeAmplitude(order, insertions, options);
}
mty::Amplitude Model::computeAmplitude(
int order,
std::vector<mty::Insertion> insertions,
FeynOptions options
Kinematics const &kinematics
)
{
Kinematics kinematics { insertions };
FeynOptions options;
return computeAmplitude(order, insertions, kinematics, options);
}
mty::Amplitude Model::computeAmplitude(
int order,
std::vector<mty::Insertion> insertions,
FeynOptions &options
)
{
Kinematics kinematics { insertions };
auto res = computeAmplitude(order, insertions, kinematics, options);
options = res.getOptions();
return res;
}
mty::Amplitude Model::computeAmplitude(
int order,
std::vector<mty::Insertion> insertions,
Kinematics const &kinematics,
FeynOptions options
FeynOptions &options
)
{
if (!options.getFeynRuleCalculation())
......@@ -108,16 +129,20 @@ mty::Amplitude Model::computeAmplitude(
if (options.getFeynRuleCalculation()) {
std::vector<Lagrangian::TermType> lagrangian = L.interaction;
options.applyFilters(lagrangian);
return computeAmplitude(
auto res = computeAmplitude(
lagrangian, insertions, kinematics, options
);
options = res.getOptions();
return res;
}
else {
std::vector<FeynmanRule> const &rules = getFeynmanRules();
auto filteredRules = options.applyFilters(rules);
return computeAmplitude(
auto res = computeAmplitude(
filteredRules, insertions, kinematics, options
);
options = res.getOptions();
return res;
}
}
......@@ -477,7 +502,7 @@ WilsonSet Model::getWilsonCoefficients(
//WilsonSet res(wilsons.begin(), wilsons.end());
//res.sort();
//res.mergeSorted();
wilsons.options = ampl.getOptions();
wilsons.options = feynOptions;
wilsons.kinematics = ampl.getKinematics();
return wilsons;
......@@ -525,8 +550,9 @@ WilsonSet Model::computeWilsonCoefficients_default(
FeynOptions const &feynOptions
)
{
auto ampl = computeAmplitude(order, insertions, feynOptions);
auto wilsons = getWilsonCoefficients(ampl, feynOptions);
auto option_cpy = feynOptions;
auto ampl = computeAmplitude(order, insertions, option_cpy);
auto wilsons = getWilsonCoefficients(ampl, option_cpy);
return wilsons;
}
......@@ -673,8 +699,10 @@ WilsonSet Model::computeWilsonCoefficients_4Fermions(
{
if (mty::option::verboseAmplitude)
std::cout << "Using special 4-fermion calculation" << std::endl;
WilsonSet res;
Kinematics kinematics { insertions };
WilsonSet res;
res.kinematics = kinematics;
res.options = feynOptions;
WilsonSet contrib;
if (feynOptions.getTopology() & Topology::Box) {
if (mty::option::verboseAmplitude)
......
#include "simpleWilson.h"
namespace mty {
}
......@@ -166,6 +166,8 @@ void WilsonSet::merge(bool sorted)
return;
}
WilsonSet other;
other.kinematics = kinematics;
other.options = options;
for (const auto &wil : *this)
addWilson(wil, other);
*this = std::move(other);
......
......@@ -17,6 +17,40 @@
#include "model.h"
#include "error.h"
namespace mty {
csl::Expr getWilsonCoefficient(
WilsonSet const &wilsons,
std::vector<Wilson> const &contributions
)
{
std::vector<csl::Expr> terms;
terms.reserve(wilsons.size());
for (const auto &wil : wilsons) {
for (const auto &contrib : contributions) {
if (wil.op == contrib.op)
terms.push_back(contrib.coef.getCoefficient()
* wil .coef.getCoefficient());
}
}
return csl::sum_s(terms);
}
csl::Expr getWilsonCoefficient(
Model const &model,
WilsonSet const &wilsons,
DiracCoupling coupling
)
{
return getWilsonCoefficient(
wilsons,
chromoMagneticOperator(model, wilsons, coupling)
);
}
} // End of namespace mty
// Old way
namespace mty {
std::ostream &operator<<(
......@@ -519,11 +553,6 @@ namespace OperatorParser {
{
HEPAssert(IsOfType<QuantumField>(qField),
mty::error::TypeError,
"This function expects a quantum field, " + toString(qField)
+ " given.")
QuantumField const *f = ConvertToPtr<QuantumField>(qField);
......@@ -592,6 +621,5 @@ namespace OperatorParser {
}
}
}
} // End of namespace mty::OperatorParser
} // End of namespace mty
} // End of namespace OperatorParser
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