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

D-term correction in the MSSM, Penguin Patch 1.0, Wilson operator standard basis

parent a3f26aeb
......@@ -3479,25 +3479,25 @@ optional<Expr> ISum::factor(bool full) const
Expr sum = copy();
Expr product = CSL_1;
for (const auto& f : factors) {
Expr expanded = Factored(sum, f.get());
if (expanded->getType() == csl::Type::Prod) {
Expr factored = Factored(sum, f.get());
if (factored->getType() == csl::Type::Prod) {
bool toDevelop = false;
for (size_t i = 0; i != expanded->size(); ++i) {
if (expanded[i]->getType() == csl::Type::Sum) {
if (expanded[i] != f) {
for (size_t i = 0; i != factored->size(); ++i) {
if (factored[i]->getType() == csl::Type::Sum) {
if (factored[i] != f) {
toDevelop = true;
sum = expanded[i];
sum = factored[i];
}
else
product = product * expanded[i];
product = product * factored[i];
}
else
product = product * expanded[i];
product = product * factored[i];
}
if (toDevelop)
return Factored(sum) * product;
else
return expanded;
return factored;
}
}
if (product != CSL_1)
......@@ -3702,7 +3702,7 @@ Expr IProd::suppressTerm(Expr_info term) const
}
}
return prod_s(newArgs);
return prod_s(newArgs, true);
}
void IProd::setArgument(const Expr& t_argument, int iArg)
......
......@@ -747,28 +747,25 @@ optional<Expr> Sum::factor(Expr_info expr, bool full) const
if (arg[i]->askTerm(expr))
toFactor.push_back(i);
else {
newAbstract = sum_s(newAbstract, arg[i]);
return std::nullopt;
}
}
else
newAbstract = sum_s(newAbstract, arg[i]);
return std::nullopt;
}
}
// If the number of arguments that can be factored is greater than 2,
// We factor them
Expr copyExpr = expr->copy();
if (toFactor.size() == argument.size()) {
Expr factored = CSL_0;
std::vector<csl::Expr> factored;
factored.reserve(toFactor.size());
for (size_t i=0; i<toFactor.size(); i++) {
Expr foo = arg[toFactor[i]]->suppressTerm(expr);
if (foo*copyExpr != arg[toFactor[i]])
return nullopt;
factored = sum_s(factored, arg[toFactor[i]]->suppressTerm(expr));
factored.push_back(arg[toFactor[i]]->suppressTerm(expr));
}
newAbstract = sum_s(newAbstract, factored * copyExpr);
return newAbstract;
return csl::prod_s(sum_s(factored), copyExpr, true);
}
else
return nullopt;
......
......@@ -135,6 +135,24 @@ class Node: public Object {
Child c = std::move(child);
addChild(c);
}
/*!
* \brief Adds a child in the Node.
* \param child New child added in \b children.
*/
virtual void addChild(Child&& child) {
children.push_back(std::move(child));
}
/*!
* \brief Adds a child in the Node.
* \param child New child added in \b children.
*/
template<class T>
inline void addChild(std::unique_ptr<T>&& child) {
Child c = std::move(child);
addChild(c);
}
/*!
* \return **children.size()**, the number of children.
......@@ -252,6 +270,35 @@ class List: public Node {
}
Node::addChild(child);
}
/*!
* \brief Adds a child in the List. Checks if the specifier is the same as
* the one of the List (each element must have the same).
* \param child New child added in \b children.
*/
void addChild(Child&& child) override {
if (child->getSpecifier() != specifier) {
std::cerr << "JSONParsingError: difference specification"
<< " \"" << child->getSpecifier() << "\" in list.\n";
exit(1);
}
Node::addChild(std::move(child));
}
/*!
* \brief Adds a child in the List. Checks if the specifier is the same as
* the one of the List (each element must have the same).
* \param child New child added in \b children.
*/
template<class T>
void addChild(std::unique_ptr<T>&& child) {
if (child->getSpecifier() != specifier) {
std::cerr << "JSONParsingError: difference specification"
<< " \"" << child->getSpecifier() << "\" in list.\n";
exit(1);
}
Node::addChild(std::move(child));
}
};
/*!
......
......@@ -135,6 +135,24 @@ class Node: public Object {
Child c = std::move(child);
addChild(c);
}
/*!
* \brief Adds a child in the Node.
* \param child New child added in \b children.
*/
virtual void addChild(Child&& child) {
children.push_back(std::move(child));
}
/*!
* \brief Adds a child in the Node.
* \param child New child added in \b children.
*/
template<class T>
inline void addChild(std::unique_ptr<T>&& child) {
Child c = std::move(child);
addChild(c);
}
/*!
* \return **children.size()**, the number of children.
......@@ -252,6 +270,35 @@ class List: public Node {
}
Node::addChild(child);
}
/*!
* \brief Adds a child in the List. Checks if the specifier is the same as
* the one of the List (each element must have the same).
* \param child New child added in \b children.
*/
void addChild(Child&& child) override {
if (child->getSpecifier() != specifier) {
std::cerr << "JSONParsingError: difference specification"
<< " \"" << child->getSpecifier() << "\" in list.\n";
exit(1);
}
Node::addChild(std::move(child));
}
/*!
* \brief Adds a child in the List. Checks if the specifier is the same as
* the one of the List (each element must have the same).
* \param child New child added in \b children.
*/
template<class T>
void addChild(std::unique_ptr<T>&& child) {
if (child->getSpecifier() != specifier) {
std::cerr << "JSONParsingError: difference specification"
<< " \"" << child->getSpecifier() << "\" in list.\n";
exit(1);
}
Node::addChild(std::move(child));
}
};
/*!
......
......@@ -14,8 +14,7 @@ namespace mty {
struct MixingInfo {
Mixing mixing;
size_t i;
size_t j;
std::vector<size_t> mixedParticles;
};
NMFV_Model(
......@@ -38,6 +37,14 @@ namespace mty {
}
void getToLowEnergyLagrangian();
void addPairMixing(
Mixing mix,
size_t i,
size_t j
);
void addMixing(MixingInfo mix);
};
......
#pragma once
#include <functional>
#include <csl.h>
namespace mty {
template<class T, class U>
class Cache {
public:
using CachedData = std::pair<T, U>;
using Comparator = std::function<bool(T const&, T const&)>;
using Releaser = std::function<U(U const&)>;
IMPLEMENTS_STD_VECTOR(CachedData, cached)
explicit
Cache(
Comparator &&t_comparator = [](T const &a, T const &b) { return a == b; },
Releaser &&t_releaser = [](U const &res) { return res; }
)
:comparator(std::forward<Comparator>(t_comparator)),
releaser(std::forward<Releaser>(t_releaser))
{
}
template<class Calculator>
U calculate(T const &input, Calculator &&calc)
{
if (auto pos = find(input) ; pos != end()) {
return releaser(pos->second);
}
return releaser(add(input, calc(input)));
}
auto find(T const &key) {
return std::find_if(cached.begin(), cached.end(),
[&](CachedData const &data) {
return comparator(data.first, key);
});
}
auto find(T const &key) const {
return std::find_if(cached.begin(), cached.end(),
[&](CachedData const &data) {
return comparator(data.first, key);
});
}
private:
U const &add(T const &key, U const &value) {
if (auto pos = find(key) ; pos == end()) {
cached.push_back({key, value});
return cached.back().second;
}
else {
return pos->second;
}
}
private:
std::vector<CachedData> cached;
Comparator comparator;
Releaser releaser;
};
} // namespace mty
......@@ -177,13 +177,6 @@ namespace mty {
*/
void setWilsonOperatorCoefficient(csl::Expr const &factor);
/**
* @brief Sets the Wilson operator basis #wilsonOperatorBasis.
*
* @param basis Basis used in Wilson coefficient decomposition.
*/
void setWilsonOperatorBasis(OperatorBasis basis);
/**
* @brief Tells if a diagram passes all diagram filters.
*
......@@ -390,6 +383,16 @@ namespace mty {
*/
void initDefaultFilters();
/**
* @brief Sets the Wilson operator basis #wilsonOperatorBasis.
*
* @note This method is private because the operator basis is for now
* fixed to the standard basis, more suitable for identification.
*
* @param basis Basis used in Wilson coefficient decomposition.
*/
void setWilsonOperatorBasis(OperatorBasis basis);
public:
/**
......
......@@ -54,6 +54,7 @@ class Lagrangian {
friend class ModelBuilder;
friend class MSSM_Model;
friend class PMSSM_Model;
friend class NMFV_Model;
friend class Expander;
public:
......
......@@ -460,6 +460,16 @@ public:
std::vector<csl::Expr> const &massMatrix
);
std::vector<std::vector<mty::Particle>> const &getParticleFamilies() const {
return particleFamilies;
}
void addParticleFamily(std::vector<mty::Particle> const &families);
void removeParticleFamily(mty::Particle const &particle);
void addParticleFamily(std::vector<std::string> const &familyNames);
void removeParticleFamily(std::string const &particleName);
protected:
void replaceTermInLagrangian(
......@@ -498,6 +508,8 @@ protected:
mty::Particle& gaugeBoson
);
void clearProjectorsInMass();
void doPromoteToMajorana(
mty::Particle &weylFermion,
std::string const &newParticleName = ""
......@@ -569,6 +581,8 @@ protected:
* mty::Library::generateSpectrum() function.
*/
std::vector<csl::Expr> abbreviatedMassExpressions;
std::vector<std::vector<mty::Particle>> particleFamilies;
};
///////////////////////////////////////////////////
......
#pragma once
#include <vector>
namespace mty {
class Wilson;
class Kinematics;
class Amplitude;
/**
* @brief Tells if a set of wilson coefficients requires the penguin patch
* implemented in this file. This function must return true for the
* applyPenguinPatch() function to work (otherwise it raises an error).
*
* @param amplitude Amplitude of the process.
*
* @return \b True if the patch is required.
* @return \b False otherwise.
*/
bool requiresPenguinPatch(
Amplitude const &amplitude
);
/**
* @brief Applies the patch for the on-shell calculation of penguins with
* massless vector.
*
* @details When the outgoing vector \f$ A \f$ of \f$ \psi\to\psi A \f$
* transition is massless, the calculation of Wilson coefficients requires
* the on-shell calculation with the inclusion of counter-terms for the
* corrections of external fermions. As these counter-terms are not taken
* care of automatically for now, this function applies gauge invariance to
* recover the correct result, using the fact that the result is redundant
* and that all quantities can be derived from the result without
* counter-terms. For box diagrams or penguins with massive vector boson,
* the calculation suffers from no issue.
*
* @param wilsons Wilson Operator / Coefficients to correct.
* @param kinematics Kinematica of the process.
*/
void applyPenguinPatch(
std::vector<Wilson> &wilsons,
Kinematics const &kinematics
);
} // namespace mty
......@@ -25,6 +25,7 @@
#include <csl.h>
#include "kinematics.h"
#include "cache.h"
#include "feynOptions.h"
namespace mty {
......@@ -136,6 +137,18 @@ struct WilsonSet: public std::vector<Wilson> {
FeynOptions options {};
};
std::vector<Wilson> copyWilsons(std::vector<Wilson> const &wilsons);
inline Cache<
csl::Expr, // Key type
std::vector<Wilson> // Result type
> cachedWilson(
[](csl::Expr const &A, csl::Expr const &B) {
return (A == B) || mty::hardComparison(A, B);
}, // Comparator
copyWilsons // Releaser
);
void parseStructures(
csl::Expr &arg,
std::vector<csl::Expr> &inOperator,
......@@ -146,6 +159,23 @@ std::vector<csl::Expr> parseStructures(
csl::Expr& product
);
std::vector<Wilson> cachedWilsonCalculation(
csl::Expr const& product,
csl::Expr const& operatorFactor,
Wilson res,
csl::Expr op,
bool standardBasis,
bool recursive
);
std::vector<Wilson> sglSimplifyForWilson(
csl::Expr const &op,
csl::Expr const &operatorFactor,
Wilson res,
bool standardBasis,
bool recursive
);
std::vector<Wilson> parseExpression(
csl::Expr const& product,
csl::Expr const& operatorFactor = CSL_1,
......
......@@ -116,7 +116,7 @@ namespace OperatorParser {
Insertion const &outgoingFermion,
Insertion const &vectorBoson,
csl::Tensor momentum,
Chirality chirality
bool chiralOperator
);
std::optional<Insertion> getInsertion(
......
......@@ -213,20 +213,6 @@ protected:
void initInteractions164();
void initInteractions165();
void initInteractions166();
void initInteractions167();
void initInteractions168();
void initInteractions169();
void initInteractions170();
void initInteractions171();
void initInteractions172();
void initInteractions173();
void initInteractions174();
void initInteractions175();
void initInteractions176();
void initInteractions177();
void initInteractions178();
void initInteractions179();
void initInteractions180();
void initSpectrum();
protected:
......
......@@ -923,20 +923,6 @@ void PMSSM_LEM::initInteractions()
initInteractions164();
initInteractions165();
initInteractions166();
initInteractions167();
initInteractions168();
initInteractions169();
initInteractions170();
initInteractions171();
initInteractions172();
initInteractions173();
initInteractions174();
initInteractions175();
initInteractions176();
initInteractions177();
initInteractions178();
initInteractions179();
initInteractions180();
}
void PMSSM_LEM::initSpectrum()
......
This diff is collapsed.
......@@ -722,7 +722,7 @@ void MSSM_Model::initQuarticDTerms()
});
csl::Replace(quadratic2, B, A);
addLagrangianTerm(
csl::Refreshed(CSL_M_HALF
csl::Refreshed(CSL_HALF
* quadratic1 * quadratic2
));
}
......@@ -1263,7 +1263,7 @@ void MSSM_Model::adjustCouplingConstants()
csl::Expr thetaW = sm_input::theta_W;
replace(
g_L*g_L + g_Y*g_Y,
csl::pow_s(g_L/csl::cos_s(thetaW), csl::int_s(2)));
csl::pow_s(g_L, 2) * (1 + csl::pow_s(csl::tan_s(thetaW), 2)));
replace(
g_Y,
e / csl::cos_s(thetaW));
......
......@@ -33,18 +33,110 @@ NMFV_Model::NMFV_Model(
checkHermiticity();
}
void NMFV_Model::getToLowEnergyLagrangian()
{
InteractionTerm::abbreviateFactors = true;
std::cout << "Breaking SU(2)_L symmetry ...\n";
breakSU2LGaugeSymmetry();
replaceWBoson();
std::cout << "Expanding around Higgses VEV's ...\n";
expandAroundVEVs();
L.mergeTerms();
std::cout << "Diagonalizing SM mass matrices ...\n";
diagonalize2By2Matrices();
diagonalizeYukawas();
L.mergeTerms();
adjustCouplingConstants();
std::cout << "Breaking flavor symmetry ...\n";
breakSMFlavorSymmetry();
approximateQuarkMasses();
approximateCKM();
std::cout << "Diagonalizing MSSM mass matrices ...\n";
L.mergeTerms();
approximateSFermionMixings();
std::cout << "Diagonalizing Neutralinos ..." << std::endl;
diagonalizeNeutralinos();
std::cout << "Diagonalizing Charginos ..." << std::endl;
diagonalizeCharginos();
std::cout << "Promoting Majorana fermions ..." << std::endl;
promoteMajoranas();
// Partial NMFV without quartic scalar interactions
for (size_t i = 0; i != L.interaction.size(); ++i) {
if (L.interaction[i]->getContent().size() == 4) {
L.interaction.erase(L.interaction.begin() + i);
--i;
}
}
std::cout << "Diagonalizing SFermions ..." << std::endl;
diagonalizeSFermions();
renameSFermions();
generateDiracFermions();
gatherMasses();
promoteToGoldstone("Gp", "W");
promoteToGoldstone("G0", "Z");
refresh();
}
void NMFV_Model::addPairMixing(
Mixing mix,
size_t i,
size_t j
)
{
using name_container = std::array<const char*, 3>;
constexpr name_container ULname = { "su_L", "sc_L", "st_L" };
constexpr name_container URname = { "su_R", "sc_R", "st_R" };
constexpr name_container DLname = { "sd_L", "ss_L", "sb_L" };
constexpr name_container DRname = { "sd_R", "ss_R", "sb_R" };