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

Divergences solved, fermion propagator corrected

parent 919f0e61
......@@ -622,6 +622,8 @@ class IndexStructure{
*/
IndexStructure operator+(const IndexStructure& structure) const;
bool hasCommonIndex(csl::IndexStructure const &other) const;
/*!
* \brief Compares the IndexStructure with \b structure. Each index must
* match exactly (see Index::exactMatch()) with the Index at the same place
......
......@@ -579,6 +579,15 @@ IndexStructure IndexStructure::operator+(const IndexStructure& structure) const
return newStructure;
}
bool IndexStructure::hasCommonIndex(const IndexStructure& structure) const
{
for (const auto &A : *this)
for (const auto &B : structure)
if (A == B)
return true;
return false;
}
bool IndexStructure::exactMatch(const IndexStructure& structure) const
{
if (index.size() != structure.size())
......
......@@ -93,6 +93,13 @@ namespace mty {
*/
int getTopology() const { return topology; }
/**
* @return #fermionOrder
*/
std::vector<int> const &getFermionOrder() const {
return fermionOrder;
}
/**
* @return #lfilters
*/
......@@ -152,6 +159,16 @@ namespace mty {
*/
void setTopology(int t_topology);
/**
* @brief Sets the fermion order used to express the amplitude.
*
* @param order Fermion order to apply, range of integers from \f$ 0 \f$
* to \f$ N-1 \f$ for \f$ N \f$ external fermions in the process.
*
* @sa #fermionOrder
*/
void setFermionOrder(std::vector<int> const &order);
/**
* @brief Sets the Wilson operator coefficient
* #wilsonOperatorCoefficient.
......@@ -456,6 +473,27 @@ namespace mty {
*/
int topology = mty::Topology::Any;
/**
* @brief Order for fermions in bilinears.
*
* @details The order of fermions in the result in fixed, useful for
* Wilson coefficient extraction. For example, a process with for
* external fermions \f$ \psi_0,\ \psi_1,\ \psi_2 \f$ and \f$ \psi_3 \f$
* can be expressed with the order \f$ (0123) \f$ with terms of the type
* \f[
* \left(\bar{\psi}_0\Gamma^A\psi_1\right)
* \left(\bar{\psi}_1\Gamma^B\psi_2\right,
* \f]
* while an order \f$ (2013) \f$ corresponds to something such as
* \f[
* \left(\bar{\psi}_2\Gamma^C\psi_0\right)
* \left(\bar{\psi}_1\Gamma^D\psi_3\right,
* \f]
* with \f$ \Gamma^{A,B,C,D} \f$ arbitrary \f$ \gamma \f$-matrix
* combinations.
*/
std::vector<int> fermionOrder;
/**
* @brief List of lagrangian filters.
*
......
......@@ -100,6 +100,8 @@ class Insertion {
csl::Expr getExpression() const;
bool isEquivalent(Insertion const &other) const;
bool operator==(Insertion const &other) const;
bool operator!=(Insertion const &other) const {
......
......@@ -194,6 +194,8 @@ public:
name = t_name;
}
bool isName(std::string_view t_name) const;
/**
* @brief Returns the first element (if found) with id \b id.
*
......@@ -445,6 +447,9 @@ std::ostream &operator<<(
LHAFileData const &data
);
std::string tolower(std::string const &str);
} // End of namespace mty::lha
#endif
......@@ -121,45 +121,15 @@ public:
FeynOptions options = {}
);
mty::Amplitude computeAmplitude(
int order,
std::vector<mty::Insertion> insertions,
std::vector<int> const &fermionOrder,
FeynOptions options = {}
);
mty::Amplitude computeAmplitude(
int order,
std::vector<mty::Insertion> insertions,
std::vector<int> const &fermionOrder,
Kinematics const &kinematics,
FeynOptions options = {}
);
mty::Amplitude computePartialAmplitude(
int order,
std::vector<mty::Insertion> insertions,
FeynOptions options = {}
);
mty::Amplitude computePartialAmplitude(
int order,
std::vector<mty::Insertion> insertions,
Kinematics const &kinematics,
FeynOptions options = {}
);
mty::Amplitude computePartialAmplitude(
int order,
std::vector<mty::Insertion> insertions,
std::vector<int> const &fermionOrder,
FeynOptions options = {}
);
mty::Amplitude computePartialAmplitude(
int order,
std::vector<mty::Insertion> insertions,
std::vector<int> const &fermionOrder,
Kinematics const &kinematics,
FeynOptions options = {}
);
......@@ -184,7 +154,6 @@ public:
mty::Amplitude computeAmplitude(
std::vector<Lagrangian::TermType> &lagrangian,
std::vector<Insertion> insertions,
std::vector<int> fermionOrder,
Kinematics const &kinematics,
FeynOptions const &options,
std::vector<FeynmanRule const*> rules = {}
......@@ -193,7 +162,6 @@ public:
mty::Amplitude computeAmplitude(
std::vector<FeynmanRule const*> &feynRules,
std::vector<Insertion> const &insertions,
std::vector<int> const &fermionOrder,
Kinematics const &kinematics,
FeynOptions const &options
);
......@@ -208,9 +176,26 @@ public:
bool applyDegreesOfFreedomFactor = true
);
csl::Expr computeSquaredAmplitude(
Amplitude const &amplL,
Amplitude const &amplR,
bool applyDegreesOfFreedomFactor = true
);
csl::Expr computeSquaredAmplitude(
WilsonSet const &amplL,
WilsonSet const &amplR,
bool applyDegreesOfFreedomFactor = true
);
WilsonSet getWilsonCoefficients(
Amplitude const &ampl,
FeynOptions const &feynOptions,
bool squaredAfter = false
);
WilsonSet getWilsonCoefficients(
Amplitude const &ampl,
FeynOptions const &feynOptions = {},
bool squaredAfter = false
);
......@@ -323,6 +308,8 @@ protected:
std::vector<mty::FeynmanRule> feynmanRules;
};
int operatorDegeneracy(std::vector<mty::Insertion> const &insertions);
} // End of namespace mty
#endif /* MODE_H_INCLUDED */
......@@ -39,7 +39,7 @@ namespace mty::option {
///////////////////////////////////////////////////
inline bool simplifyAmplitudes = true;
inline bool orderExternalFermions = false;
inline bool orderExternalFermions = true;
inline bool discardLowerOrders = true;
inline bool evaluateFermionTraces = true;
inline bool excludeTadpoles = true;
......
......@@ -13,98 +13,12 @@
// You should have received a copy of the GNU General Public License
// 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';
}
#include <marty>
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)
......@@ -115,6 +29,25 @@ int main() {
mty::TwoHDM_Model<type> THDM;
std::cout << THDM << std::endl;
FeynOptions options;
options.setTopology(Topology::Box);
options.addFilter(filter::forceParticleCombination({"W", "t"}));
options.addFilter(filter::disableParticles({"u_L", "c_L", "H^+", "G^+"}));
options.setFermionOrder({3, 0, 1, 2});
mty::option::keepOnlyFirstMassInLoop = true;
auto ampl = THDM.computeAmplitude(
OneLoop,
{Incoming("d"), Incoming(AntiPart("s")), Outgoing(AntiPart("d")), Outgoing("s")},
options
);
Display(ampl);
Show(ampl);
auto wilsons = THDM.getWilsonCoefficients(ampl);
for (auto &coef : wilsons)
coef.coef.setCoefficient(csl::Evaluated(coef.coef.getCoefficient(), csl::eval::abbreviation));
Display(wilsons);
return 0;
}
......@@ -134,10 +134,8 @@ namespace mty {
)
{
std::vector<mty::Insertion*> fermions = fermionsOf(insertions);
HEPAssert(fermions.size() == order.size(),
mty::error::ValueError,
"Found " + toString(fermions.size()) + " for a fermion ordering"
" of " + toString(order.size()) + " fermions.")
if (fermions.size() != order.size())
return;
for (size_t i = 0; i != order.size()/2; ++i) {
const int i1 = order[2*i];
const int i2 = order[2*i + 1];
......
......@@ -37,6 +37,11 @@ namespace mty {
topology = t_topology;
}
void FeynOptions::setFermionOrder(std::vector<int> const &order)
{
fermionOrder = order;
}
void FeynOptions::setWilsonOperatorCoefficient(csl::Expr const &factor)
{
wilsonOperatorCoefficient = csl::DeepCopy(factor);
......
......@@ -111,6 +111,14 @@ csl::Expr Insertion::getExpression() const
return insertion;
}
bool Insertion::isEquivalent(Insertion const &other) const
{
if (field != other.field || mediator != other.mediator)
return false;
return field->isSelfConjugate()
|| ((particle == other.particle) == (incoming == other.incoming));
}
bool Insertion::operator==(Insertion const &other) const
{
return field == other.field
......
......@@ -806,6 +806,45 @@ ostream& operator<<(ostream & out,
return out;
}
static void sortTensors(std::vector<csl::Expr> &tensors)
{
auto free = [&](csl::IndexStructure const &index) {
return csl::Abbrev::getFreeStructure(index);
};
std::sort(tensors.begin(), tensors.end());
std::reverse(tensors.begin(), tensors.end());
std::vector<csl::Expr> sorted;
sorted.reserve(tensors.size());
csl::IndexStructure contractedIndices;
auto step = [&](size_t pos) {
contractedIndices = free(
contractedIndices + tensors[pos]->getIndexStructureView());
sorted.push_back(tensors[pos]);
tensors.erase(tensors.begin() + pos);
};
while (!tensors.empty()) {
if (contractedIndices.empty()) {
step(0);
continue;
}
bool foundCommon = false;
for (size_t i = 0; i != tensors.size(); ++i) {
if (dynamic_cast<mty::QuantumField const*>(tensors[i].get()))
continue;
auto const &index = tensors[i]->getIndexStructureView();
if (contractedIndices.hasCommonIndex(index)) {
step(i);
foundCommon = true;
break;
}
}
if (!foundCommon)
step(0);
}
tensors = std::move(sorted);
}
int matchBOnA(csl::Expr const& A, csl::Expr &B)
{
std::vector<csl::Expr> tensorsInA;
......@@ -823,10 +862,10 @@ int matchBOnA(csl::Expr const& A, csl::Expr &B)
if (tensorsInA.size() != tensorsInB.size()) {
return tensorsInA.size() < tensorsInB.size();
}
std::sort(tensorsInA.begin(), tensorsInA.end());
std::sort(tensorsInB.begin(), tensorsInB.end());
sortTensors(tensorsInA);
sortTensors(tensorsInB);
std::vector<std::pair<csl::Index, csl::Index>> mapping;
for (size_t i = tensorsInA.size(); i --> 0 ;)
for (size_t i = 0; i != tensorsInA.size(); ++i) {
if (tensorsInA[i]->getParent_info()
!= tensorsInB[i]->getParent_info()) {
return tensorsInA[i]->getName() < tensorsInB[i]->getName();
......@@ -852,7 +891,7 @@ int matchBOnA(csl::Expr const& A, csl::Expr &B)
mapping.push_back({ Bstruct[j], Astruct[j] });
}
}
}
std::vector<csl::Index> intermediateIndices;
intermediateIndices.reserve(mapping.size());
for (const auto &mappy : mapping)
......
......@@ -87,22 +87,22 @@ namespace mty {
csl::Expr Kinematics::getDegreesOfFreedomFactor() const
{
csl::Expr factor = CSL_1;
std::vector<size_t> incomingIndices;
incomingIndices.reserve(insertions.size());
std::vector<size_t> outgoingIndices;
outgoingIndices.reserve(insertions.size());
for (size_t i = 0; i != insertions.size(); ++i) {
if (insertions[i].isIncoming())
factor *= insertions[i].getField()->getNDegreesOfFreedom();
else
incomingIndices.push_back(i);
outgoingIndices.push_back(i);
}
while (!incomingIndices.empty()) {
const size_t i = incomingIndices[0];
while (!outgoingIndices.empty()) {
const size_t i = outgoingIndices[0];
std::vector<size_t> identicalFields;
identicalFields.reserve(insertions.size());
identicalFields.push_back(0);
for (size_t index = 1; index != incomingIndices.size(); ++index) {
const size_t j = incomingIndices[index];
for (size_t index = 1; index != outgoingIndices.size(); ++index) {
const size_t j = outgoingIndices[index];
if (areIdenticalFields(insertions[i], insertions[j])) {
identicalFields.push_back(index);
}
......@@ -110,7 +110,7 @@ namespace mty {
int combi = std::tgamma(static_cast<int>(1+identicalFields.size()));
if (combi != 1)
factor *= combi;
removeIndices(incomingIndices, identicalFields);
removeIndices(outgoingIndices, identicalFields);
}
return factor;
}
......
......@@ -17,6 +17,7 @@
#include "lhaData.h"
#include <iostream>
#include <fstream>
#include <algorithm>
namespace mty::lha {
......@@ -105,7 +106,7 @@ std::vector<std::string> Reader::split(std::string const &line)
bool Reader::isBlockLine(std::vector<std::string> const &line)
{
return line.size() >= 2
and line[0] == "BLOCK";
and mty::lha::tolower(line[0]) == "block";
}
bool Reader::isPositiveInteger(std::string const &str)
......
......@@ -91,6 +91,17 @@ LHABlock::LHABlock(std::string const &t_name)
setName(t_name);
}
bool LHABlock::isName(std::string_view t_name) const
{
if (name.size() != t_name.size())
return false;
for (size_t i = 0; i != name.size(); ++i) {
if (std::tolower(name[i]) != std::tolower(t_name[i]))
return false;
}
return true;
}
std::optional<LHAElement> LHABlock::getElement(size_t id) const
{
for (const auto &e : elements) {
......@@ -174,7 +185,7 @@ size_t LHAFileData::findBlock(
) const
{
for (size_t i = 0; i != blocks.size(); ++i)
if (blocks[i].getName() == nameBlock)
if (blocks[i].isName(nameBlock))
return i;
return npos;
}
......@@ -191,7 +202,7 @@ std::vector<FloatType> LHAFileData::getValues(
) const
{
for (const auto &block : blocks)
if (block.getName() == nameBlock) {
if (block.isName(nameBlock)) {
std::vector<FloatType> res;
res.reserve(block.size());
for (const auto &el : block) {
......@@ -209,7 +220,7 @@ std::optional<FloatType> LHAFileData::getValue(
) const
{
for (const auto &block : blocks)
if (block.getName() == nameBlock)
if (block.isName(nameBlock))
for (const auto &el : block)
if (el.id == id)
return el.value;
......@@ -224,7 +235,7 @@ std::optional<FloatType> LHAFileData::getValue(
) const
{
for (const auto &block : blocks)
if (block.getName() == nameBlock)
if (block.isName(nameBlock))
for (const auto &el : block)
if (el.id == i and el.id_sup == j)
return el.value;
......@@ -290,4 +301,13 @@ std::ostream &operator<<(
return out;
}
std::string tolower(std::string const &str)
{
std::string cpy(str);
std::for_each(cpy.begin(), cpy.end(), [&](char &c) {
c = std::tolower(c);
});
return cpy;
}
} // End of namespace lha
......@@ -102,60 +102,20 @@ mty::Amplitude Model::computeAmplitude(
FeynOptions options
)
{
auto fermionOrder = defaultFermionOrder(insertions);
if (!options.getFeynRuleCalculation())
options.setLoopOrder(order, static_cast<int>(insertions.size()));
if (options.getFeynRuleCalculation()) {
std::vector<Lagrangian::TermType> lagrangian = L.interaction;
options.applyFilters(lagrangian);
return computeAmplitude(
lagrangian, insertions, fermionOrder, kinematics, options
lagrangian, insertions, kinematics, options
);
}
else {
std::vector<FeynmanRule> const &rules = getFeynmanRules();
auto filteredRules = options.applyFilters(rules);
return computeAmplitude(
filteredRules, insertions, fermionOrder, kinematics, options
);
}
}
mty::Amplitude Model::computeAmplitude(