Updated MacOS installation

.setup_script.sh

@@ -126,7 +126,12 @@ testDependencies_Linux()
 testDependency_MACOS()
 {
     echo "Searching for $1 ..."
-    brew list --$2 $1 &>/dev/null || brew $2 install $1
+    if [ $# -ge 2 ]
+    then
+        brew list --$2 $1 &>/dev/null || brew install --$2 $1
+    else
+        brew list $1 &>/dev/null || brew install $1
+    fi
     brew cleanup
 }
 

csl/include/librarygroup.h

@@ -113,6 +113,18 @@ namespace csl {
                 int nIndent
                 ) const;
 
+        void printNameMapElement(
+                std::ostream       &out,
+                int                 nIndent,
+                LibParameter const &param,
+                std::string  const &type
+                ) const;
+
+        void printNameMap(
+                std::ostream &out,
+                int           nIndent
+                ) const;
+
         void printStructDefinition(
                 std::ostream &out,
                 int nIndent

csl/src/abreviation.cpp

@@ -503,7 +503,8 @@ Expr Abbrev::makeAbbreviation(std::string name,
     if (encaps->size() == 0) // nothing to abbreviate
         return encaps;
     if (csl::IsProd(encaps) && csl::IsNumerical(encaps[0])
-            && csl::Size(encaps) > 2) {
+            && csl::Size(encaps) > 2
+            && getAbbreviationsForName(name).size() > 100) {
         auto prod = csl::prod_s(
                 std::vector<csl::Expr>(encaps->begin()+1, encaps->end()), true);
         auto prodAbbrev = makeAbbreviation(name, prod, split);

csl/src/hardComparison.cpp

@@ -134,7 +134,8 @@ static bool hardComparison_impl(
     const int match = matchBOnA(A, B);
     if (match != -1)
         return false;
-    return A->compareWithDummy(B.get());
+    auto res = A->compareWithDummy(B.get());
+    return res;
 }
 
 bool hardComparison(

csl/src/librarygenerator.cpp

@@ -288,6 +288,7 @@ namespace csl {
             file paramFile(path + "/" + incDir + "/params.h");
             paramFile << "#ifndef CSL_LIB_PARAM_H_INCLUDED\n",
             paramFile << "#define CSL_LIB_PARAM_H_INCLUDED\n\n";
+            paramFile << "#include <map>\n";
             paramFile << "#include <array>\n";
             paramFile << "#include \"common.h\"\n";
             paramFile << "#include \"libcomplexop.h\"\n";

csl/src/librarygroup.cpp

@@ -270,9 +270,55 @@ namespace csl {
         printResetDefinition(out, nIndent + 1);
         out << '\n';
         printPrintDefinition(out, nIndent + 1);
+        out << '\n';
+        printNameMap(out, nIndent + 1);
         out << LibraryGenerator::indent(nIndent) << "};\n";
     }
 
+    void LibraryGroup::printNameMapElement(
+            std::ostream       &out,
+            int                 nIndent,
+            LibParameter const &param,
+            std::string  const &type
+            ) const
+    {
+        if (param.type == type) {
+            out << LibraryGenerator::indent(1+nIndent) 
+                << "{\"" << param.name << "\", &" << param.name << "},\n";
+        }
+    }
+
+    void LibraryGroup::printNameMap(
+            std::ostream &out,
+            int           nIndent
+            ) const
+    {
+        out << LibraryGenerator::indent(nIndent)
+            << "std::map<std::string, csl::InitSanitizer<"
+            << LibraryGenerator::realUsing << ">*> realParams {\n";
+        for (const auto &param : parameters) {
+            printNameMapElement(
+                    out, nIndent, param, LibraryGenerator::realUsing);
+        }
+        for (const auto &param : forcedParameters) {
+            printNameMapElement(
+                    out, nIndent, param, LibraryGenerator::realUsing);
+        }
+        out << "};\n\n";
+        out << LibraryGenerator::indent(nIndent)
+            << "std::map<std::string, csl::InitSanitizer<" 
+            << LibraryGenerator::complexUsing << ">*> complexParams {\n";
+        for (const auto &param : parameters) {
+            printNameMapElement(
+                    out, nIndent, param, LibraryGenerator::complexUsing);
+        }
+        for (const auto &param : forcedParameters) {
+            printNameMapElement(
+                    out, nIndent, param, LibraryGenerator::complexUsing);
+        }
+        out << LibraryGenerator::indent(nIndent) << "};\n\n";
+    }
+
     void LibraryGroup::printFunctionStack(
             std::ostream &out,
             int           nIndent
@@ -288,6 +334,18 @@ namespace csl {
             out << "Callable{\"" << f.getName() << "\", " << f.getName() << "},\n";
         }
         out << LibraryGenerator::indent(nIndent) << "};\n\n";
+
+        out << LibraryGenerator::indent(nIndent)
+            << "inline std::map<std::string, Callable<" << (complexReturn ? 
+                    LibraryGenerator::complexUsing 
+                    :LibraryGenerator::realUsing) << ", " << getParamName();
+        out << ">> fmap_" << name << " {\n";
+        for (size_t i = 0; i != functions.size(); ++i) {
+            out << LibraryGenerator::indent(nIndent + 1);
+            out << "{\"" << functions[i].getName() 
+                << "\", f_" << name << "[" << i << "]},\n";
+        }
+        out << LibraryGenerator::indent(nIndent) << "};\n\n";
     }
 
     void LibraryGroup::printParameterDefinition(

grafed/grafed-gui/diagramrenderer.cpp

@@ -32,6 +32,8 @@ DiagramRenderer::DiagramRenderer(QWidget *parent) :
     nColumns(4),
     minWidth(300),
     minHeight(300),
+    firstDiag(0),
+    lastDiag(0),
     edgeForm(nullptr),
     nodeForm(nullptr),
     zoomValue(1.),
@@ -303,6 +305,24 @@ void DiagramRenderer::resetSelection()
     emit nonEmptySelection(false);
 }
 
+// Function not ready to be used
+void DiagramRenderer::setDiagramSet(qint32 first, qint32 last)
+{
+    if (mode == SessionMode) {
+        firstDiag = first;
+        lastDiag  = last;
+        diagrams = QList<Diagram*>{};
+        table->clear();
+        for (qint32 i = first; i != last; ++i) {
+            auto res = newDiagram();
+            res.first->loadLinker(links[static_cast<size_t>(i)]);
+            res.first->showNodes(false);
+            diagrams.append(res.first);
+        }
+        init();
+    }
+}
+
 void DiagramRenderer::closeNodeForm()
 {
     if (nodeForm) {
@@ -408,6 +428,8 @@ std::pair<Diagram*, DiagramWidget*> DiagramRenderer::newDiagram()
     }
     auto [diagram, graph] = generateDiagram(i, j);
     diagrams.push_back(diagram);
+    allDiagrams.push_back(diagram);
+    ++lastDiag;
     table->setCellWidget(i, j, graph);
     table->cellWidget(i, j)->setToolTip(
                 QString("%1").arg(i*nColumns + j));
@@ -419,7 +441,11 @@ void DiagramRenderer::readFile(QString const &nameFile)
 {
     links = drawer::LatexLinker::loadMultiple(nameFile.toStdString());
     table->clear();
+    constexpr size_t maxDiagrams = 200;
+    size_t nDiagrams = 0;
     for (auto& link : links) {
+        if (++nDiagrams > maxDiagrams)
+            break;
         auto res = newDiagram();
         res.first->loadLinker(link);
         res.first->showNodes(false);

grafed/grafed-gui/diagramrenderer.h

@@ -48,6 +48,10 @@ class DiagramRenderer : public QWidget
         SessionMode
     };
 
+public:
+
+    static inline qint32 maxNDiagrams = 100;
+
 private:
     explicit DiagramRenderer(QWidget *parent = nullptr);
 public:
@@ -80,6 +84,9 @@ public slots:
     void cellPressed(int i, int j);
     void resetSelection();
 
+    // Function not ready to be used
+    void setDiagramSet(qint32 first, qint32 last);
+
     bool hasEdgeForm() {
         return edgeForm;
     }
@@ -134,7 +141,10 @@ private:
     qint32          nColumns;
     qint32          minWidth;
     qint32          minHeight;
+    qint32          firstDiag;
+    qint32          lastDiag;
     std::vector<drawer::LatexLinker> links;
+    QList<Diagram*>       allDiagrams;
     QList<Diagram*>       diagrams;
     QList<std::pair<int, int>> selectedDiagrams;
     qint32                indexDiagram;

grafed/grafed-gui/latexcompiler.cpp

@@ -55,7 +55,10 @@ QPixmap latexcompiler::generateLabel(const std::string &texCode)
     std::ofstream f(grafed_dir + "/eq.tex");
     if (!f)
         return QPixmap();
-    f << "\\documentclass[preview]{standalone}";
+    f << "\\documentclass[preview]{standalone}\n";
+    f << "\\usepackage{amsmath}\n";
+    f << "\\usepackage{slashed}\n";
+    f << "\\usepackage{xcolor}\n";
     f << "\\begin{document}$" << texCode << "$\\end{document}";
     f.close();
     res = system(("cd " + grafed_dir + " && "

marty/include/SM.h

@@ -36,7 +36,18 @@ namespace mty {
     class SM_Model: public mty::Model {
 
     public:
-        SM_Model();
+        SM_Model(bool init = true);
+
+        void init();
+
+        void HiggsVEVExpansion();
+        void diagonalizeSMMassMatrices();
+        void gaugeSymmetryBreaking();
+        void replaceLeptonYukawa();
+        void replaceUpYukawa();
+        void replaceDownYukawa();
+        void flavorSymmetryBreaking();
+        void adjust();
     };
 
 } // End of namespace mty

marty/include/fermionicField.h

@@ -129,6 +129,16 @@ class WeylFermion: public QuantumFieldParent {
 
     void setDiracParent(Particle const &parent) override;
 
+    void breakParticle(
+            mty::Group                     *brokenGroup,
+            std::vector<std::string> const &newNames
+            ) override;
+    void breakParticle(
+            mty::FlavorGroup                     *brokenFlavor,
+            std::vector<mty::FlavorGroup*> const &subGroups,
+            std::vector<std::string>       const &names
+            ) override;
+
     protected:
 
     void initPropagator() override;
@@ -235,9 +245,21 @@ class DiracFermion: public QuantumFieldParent {
 
     Particle getWeylFermion(Chirality chirality) const override;
 
+    void breakParticle(
+            mty::Group                     *brokenGroup,
+            std::vector<std::string> const &newNames
+            ) override;
+    void breakParticle(
+            mty::FlavorGroup                     *brokenFlavor,
+            std::vector<mty::FlavorGroup*> const &subGroups,
+            std::vector<std::string>       const &names
+            ) override;
+
     private:
 
     void initPropagator() override;
+
+    void updateChiralBrokenParts(csl::Space const *space);
 };
 
 /**

marty/include/feynOptions.h

@@ -542,6 +542,8 @@ namespace mty {
          * the standard one uses \f$ \1,\gamma^5 \f$.
          */
         OperatorBasis wilsonOperatorBasis = OperatorBasis::Chiral;
+
+        constexpr static size_t nDefaultDiagramFilters = 1;
     };
 
 } // namespace mty

marty/include/looptools_interface.h

@@ -4,6 +4,7 @@
 
 namespace mty::lt {
 
+inline bool       rescaleIntegral = false;
 inline lti_real_t minScale = 1e-3;
 inline lti_real_t regScale = 1e+0;
 
@@ -48,6 +49,7 @@ inline lti_complex_t C0iC(
         lti_real_t const &eps = 0.
         )
 {
+    if (rescaleIntegral) {
        auto min = 10*minScale; // Ensure initial value higher than minscale
        for (const auto &value: {p1, p2, p3, m1, m2, m3}) {
            const auto re = abs_re(value);
@@ -60,6 +62,7 @@ inline lti_complex_t C0iC(
                    scale, scaleFactorC(id, scale), lt::C0iC, 
                    id, eps, p1, p2, p3, m1, m2, m3);
        }
+    }
     regulateC(p1, p2, p3, m1, m2, m3, eps);
     const std::complex<lti_real_t> res = ::C0iC(
             id,

marty/include/looptools_patch.h

@@ -115,30 +115,20 @@ namespace mty::lt {
         const bool p1Zero = (p1_2 == 0.);
         const bool p3Zero = (p3_2 == 0.);
         if (p1Zero && p3Zero && p2_2 == 0.) {
-            auto min_m = abs_min_ct(m1_2, m2_2, m3_2);
-            p1_2 = eps * (min_m == 0. ? 1. : min_m);
+            p1_2 = 1e-5;
             return;
         }
-        else if (p1_2 == p2_2) {
-            if (p1Zero) {
-                p2_2 = -eps * p3_2; // p3 is different from zero here
-            }
+        else if (p1_2 == p2_2 && !p1Zero) {
             p2_2 = (1 + eps) * p1_2;
             if (abs_re(m1_2/p1_2) > 1e5 || abs_re(m2_2/p1_2) > 1e5 || abs_re(m3_2/p1_2) > 1e5)
                 p1_2 *= -1;
         }
-        else if (p1_2 == p3_2) {
-            if (p1Zero) {
-                p3_2 = -eps * p2_2; // p2 is different from zero here
-            }
+        else if (p1_2 == p3_2 && !p1Zero) {
             p3_2 = (1 + eps) * p1_2;
             if (abs_re(m1_2/p2_2) > 1e5 || abs_re(m2_2/p2_2) > 1e5 || abs_re(m3_2/p2_2) > 1e5)
                 p2_2 *= -1;
         }
-        else if (p3_2 == p2_2) {
-            if (p3Zero) {
-                p3_2 = -eps * p1_2; // p1 is different from zero here
-            }
+        else if (p3_2 == p2_2 && !p3Zero) {
             p3_2 = (1 + eps) * p2_2;
             if (abs_re(m1_2/p3_2) > 1e5 || abs_re(m2_2/p3_2) > 1e5 || abs_re(m3_2/p3_2) > 1e5)
                 p3_2 *= -1;

marty/include/mrtInterface.h

@@ -515,7 +515,6 @@ void ContractIndices(std::vector<csl::Expr> &init);
 void Display(std::vector<FeynmanRule> const& rules,
              std::ostream                  & out = std::cout);
 
-void Show(std::vector<FeynmanRule> const& rules);
 
 void Display(std::vector<csl::Expr> const& amplitudes,
              std::ostream           & out = std::cout,
@@ -527,9 +526,26 @@ void Display(mty::Amplitude const& amplitudes,
 void Display(WilsonSet const& wilsons,
              std::ostream             & out = std::cout);
 
+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(
+        std::vector<FeynmanRule> const& rules,
+        size_t first,
+        size_t last
+        );
+void Show(
+        std::vector<std::shared_ptr<wick::Graph>> const& diagrams,
+        size_t first,
+        size_t last
+        );
+void Show(
+        mty::Amplitude const& diagrams,
+        size_t first,
+        size_t last
+        );
+
 template<class T>
 void Display(T const& printable,
              std::ostream& out = std::cout)

marty/include/quantumField.h

@@ -1290,11 +1290,11 @@ class QuantumFieldParent: public csl::TensorFieldParent {
 
     operator csl::Expr();
 
-    void breakParticle(
+    virtual void breakParticle(
             mty::Group                     *brokenGroup,
             std::vector<std::string> const &newNames
             );
-    void breakParticle(
+    virtual void breakParticle(
             mty::FlavorGroup                     *brokenFlavor,
             std::vector<mty::FlavorGroup*> const &subGroups,
             std::vector<std::string>       const &names

marty/include/vectorField.h

@@ -298,9 +298,22 @@ class VectorBoson: public BaseVectorBoson {
     csl::Expr operator()(std::vector<csl::Index> indices,
                     const csl::Tensor&     t_point) override;
 
+    void breakParticle(
+            mty::Group                     *brokenGroup,
+            std::vector<std::string> const &newNames
+            ) override;
+    void breakParticle(
+            mty::FlavorGroup                     *brokenFlavor,
+            std::vector<mty::FlavorGroup*> const &subGroups,
+            std::vector<std::string>       const &names
+            ) override;
+
+
     protected:
 
     void initPropagator() override;
+
+    void updateBrokenFieldStrength(csl::Space const *space);
 };
 
 class GoldstoneBoson;

marty/src/NMFV.cpp

@@ -35,7 +35,6 @@ NMFV_Model::NMFV_Model(
 
 void NMFV_Model::getToLowEnergyLagrangian()
 {
-    InteractionTerm::abbreviateFactors = true;
     std::cout << "Breaking SU(2)_L symmetry ...\n";
     breakSU2LGaugeSymmetry();
     replaceWBoson();

marty/src/SM.cpp

@@ -26,7 +26,7 @@ using namespace csl;
 namespace mty {
 
 
-SM_Model::SM_Model()
+SM_Model::SM_Model(bool initialize)
     :mty::Model("models/files/SM.json")
 {
     getParticle("G")->setDrawType(drawer::ParticleType::Gluon);
@@ -37,9 +37,10 @@ SM_Model::SM_Model()
     csl::Index i = GaugeIndex(*this, "SU2L", H);
     csl::Index j = GaugeIndex(*this, "SU2L", H);
 
+    csl::Expr mh = sm_input::m_h;
     csl::Expr H2 = csl::GetComplexConjugate(H(i)) * H(i);
-    csl::Expr m2 = csl::constant_s("m2");
-    csl::Expr lam = csl::constant_s("lambda");
+    csl::Expr m2 = mh*mh / 2;
+    csl::Expr lam = mh*mh / (2*v*v);
 
     // Mexican hat potential
     addLagrangianTerm(m2*H2);
@@ -48,26 +49,43 @@ SM_Model::SM_Model()
     //           lam = mh^2 / (2*v^2)
     //           (With H0 -> (v + h0) / sqrt(2))
 
+    if (initialize) {
+        init();
+        refresh();
+    }
+}
+
+void SM_Model::init()
+{
+    gaugeSymmetryBreaking();
+    HiggsVEVExpansion();
+    diagonalizeSMMassMatrices();
+    replaceLeptonYukawa();
+    replaceUpYukawa();
+    replaceDownYukawa();
+    flavorSymmetryBreaking();
+    adjust();
+}
+
+void SM_Model::gaugeSymmetryBreaking()
+{
     ///////////////////////////////////////////////////
     // Breaking gauge SU(2)_L symmetry, renaming
     ///////////////////////////////////////////////////
 
     BreakGaugeSymmetry(*this, "U1Y");
-    BreakGaugeSymmetry(
-            *this,
-            "SU2L",
-            {"H", "W", "Q", "L"},
-            {{"H0", "H1"},
-            {"W1", "W2", "W3"},
-            {"U_L", "D_L"},
-            {"Nu_L", "E_L"}});
+    BreakGaugeSymmetry( *this, "SU2L");
+    renameParticle("Q_1", "U_L");
+    renameParticle("Q_2", "D_L");
+    renameParticle("L_1", "Nu_L");
+    renameParticle("L_2", "E_L");
 
     ///////////////////////////////////////////////////
     // Replacements to get SM particles W +-
     ///////////////////////////////////////////////////
 
-    Particle W1   = GetParticle(*this, "W1");
-    Particle W2   = GetParticle(*this, "W2");
+    Particle W1   = GetParticle(*this, "W_1");
+    Particle W2   = GetParticle(*this, "W_2");
     Particle W_SM = GenerateSimilarParticle("W", W1);
     W_SM->setSelfConjugate(false);
 
@@ -90,15 +108,18 @@ SM_Model::SM_Model()
     replace(
             GetFieldStrength(W2),
             CSL_I * (F_W_p - F_W_m) / csl::sqrt_s(2));
+}
 
+void SM_Model::HiggsVEVExpansion()
+{
     ///////////////////////////////////////////////////
     // Actual gauge (spontaneous) symmetry breaking
     ///////////////////////////////////////////////////
 
     csl::Expr v = sm_input::v;
 
-    Particle H0 = getParticle("H0");
-    Particle H1 = getParticle("H1");
+    Particle H1 = getParticle("H_1");
+    Particle H2 = getParticle("H_2");
 
     Particle h0 = scalarboson_s("h0 ; h^0", *this); // SM Higgs boson
     Particle Gp = scalarboson_s("Gp ; G^+", *this);
@@ -106,20 +127,19 @@ SM_Model::SM_Model()
     h0->setSelfConjugate(true);
     G0->setSelfConjugate(true);
 
-    replace(H0, Gp());
-    replace(H1, (h0() + CSL_I*G0() + v)/csl::sqrt_s(2));
-
-    csl::Expr mh = sm_input::m_h;
-    replace(m2,  mh*mh / 2);
-    replace(lam, mh*mh / (2*v*v));
+    replace(H1, Gp());
+    replace(H2, (h0() + CSL_I*G0() + v)/csl::sqrt_s(2));
+}
 
+void SM_Model::diagonalizeSMMassMatrices()
+{
     ///////////////////////////////////////////////////
     // Diagonalizing what can be
     ///////////////////////////////////////////////////
 
     DiagonalizeMassMatrices(*this);
     renameParticle("B", "A");
-    renameParticle("W3", "Z");
+    renameParticle("W_3", "Z");
 
     csl::Expr gY = getScalarCoupling("g_Y");
     csl::Expr gL = getScalarCoupling("g_L");
@@ -135,26 +155,24 @@ SM_Model::SM_Model()
     replace(
             gL, 
             e / csl::sin_s(theta_Weinberg));
+}
 
+void SM_Model::replaceLeptonYukawa()
+{
     ///////////////////////////////////////////////////
-    // Taking care of yukawa couplings
+    // Taking care of Yukawa couplings
     ///////////////////////////////////////////////////
 
     csl::Tensor Ye = GetYukawa(*this, "Ye");
-    csl::Tensor Yu = GetYukawa(*this, "Yu");
-    csl::Tensor Yd = GetYukawa(*this, "Yd");
 
     csl::Expr m_e   = sm_input::m_e;
     csl::Expr m_mu  = sm_input::m_mu;
     csl::Expr m_tau = sm_input::m_tau;
-    csl::Expr m_u   = sm_input::m_u;
-    csl::Expr m_c   = sm_input::m_c;
-    csl::Expr m_t   = sm_input::m_t;
-    csl::Expr m_d   = sm_input::m_d;
-    csl::Expr m_s   = sm_input::m_s;
-    csl::Expr m_b   = sm_input::m_b;
 
     const csl::Space* flavorSpace = GetSpace(Ye);
+    csl::Index f_i = GetIndex(flavorSpace);
+    csl::Index f_j = GetIndex(flavorSpace);