Include MARTY in C++ programs
For the samples presented here to work, one must have downloaded and
installed MARTY, for example following
instructions given on the download page.
We must here include the main header file of marty. All the library will be included in it,
including the CSL library. using namespace
instructions are not necessary, but will prevent us to write
std:: in front of objects from the C++ standard library,
mty:: in from of those in MARTY,
and csl:: for CSL.
// My first program with MARTY
#include <marty.h>
using namespace std;
using namespace mty;
using namespace csl;
int main () {
cout << "Hello World!" << endl;
return 0;
}
Compile a MARTY program
Let's now compile the code written in main.cpp:
g++ -std=c++17 -c main.cpp -o main.o
In case MARTY is installed in a non-standard path, the compiler will need to know where to find the header file (marty.h) except proper environment variables are already set (see installation instructions):
g++ -std=c++17 -c main.cpp -o main.o -I<marty-install-path>/include
One must then link the compiled file to the MARTY libraries in order to get the final executable. This gives:
g++ -std=c++17 -o main.x main.o -lmarty
In case MARTY and/or CSL are installed in a non-standard path, the compiler needs to know where to find the library files (.so files) except if proper environment variables are already set (see installation instructions):
g++ -std=c++17 -o main.x main.o -L<marty-install-path>/lib -lmarty
We are now ready to execute the compiled program:
./main.x
If MARTY and CSL are not installed in a standard path, one needs to tell the dynamic linker where to find the libraries by specifying LD_LIBRARY_PATH except if proper environment variables are already set (see installation instructions)
On Linux:
export LD_LIBRARY_PATH=<marty-install-path>/lib && ./main.x
On Mac:
export DYLD_LIBRARY_PATH=<marty-install-path>/lib && ./main.x
A first toy model example
One may download the file corresponding to this example by clicking here.
In the main.cpp, let's create with a toy model SU(2)LxU(1)Y that breaks into a U(1)em gauge. Let's consider only four Weyl fermions, QL, UR, and DR. The irreducible representations are:
- QL: (2, 1/6) SU(2) doublet with hypercharge 1/6
- UR: (1, 2/3) SU(2) singlet with hypercharge 2/3
- DR: (1, -1/3) SU(2) singlet with hypercharge -1/3
Model building
First, we have to create the model and the gauge:
Model toyModel;
toyModel.addGaugedGroup(group::Type::SU, "L", 2);
toyModel.addGaugedGroup(group::Type::U1, "Y");
toyModel.init();
We simply add the two gauged group one after the other. Once the gauge is complete,
the model must be initialized with init() before adding matter content.
Let's now name our vector bosons to follow conventions we like taking into account
the fact that MARTY set default boson names A_<group-name>:
toyModel.renameParticle("A_L", "W");
toyModel.renameParticle("A_Y", "B");
Now the matter content, the four Weyl fermions:
Particle Q_L = weylfermion_s("Q_L", toyModel, Chirality::Left);
Particle u_R = weylfermion_s("u_R", toyModel, Chirality::Right);
Particle d_R = weylfermion_s("d_R", toyModel, Chirality::Right);
Q_L->setGroupRep("L", {1});
Q_L->setGroupRep("Y", {1, 6});
u_R->setGroupRep("Y", {2, 3});
d_R->setGroupRep("Y", {-1, 3});
Several things to note here. First, one can see that when creating the fermions, one must give the model and the chirality (not for a Dirac fermion). Then, the representations must be given for each group (if not trivial). This works with the fractional charge for U(1), and dinkyn labels for the other groups. Numerator and denominator of the U(1) charge are given separately in braces, to match the syntax for dinkyn labels also given between braces (only one label here for SU(2)). For Su(2), the doublet representation has label 1, the triplet has label 2 etc. Once the representation of a particle is complete, it may be added to the model to initialize its kinetic (and possibly mass) term automatically.
toyModel.addParticle(Q_L);
toyModel.addParticle(u_R);
toyModel.addParticle(d_R);
We can now display the model in the terminal using standard output:
cout << toyModel << endl;
Model breaking
Let's now break the symmetry. For U(1)Y then SU(2)L:
toyModel.breakGaugeSymmetry("Y");
toyModel.breakGaugeSymmetry("L");
Breaking the hypercharge symmetry does basically nothing but telling MARTY that future interaction terms may violate it. Breaking the SU(2) symmetry does more. It breaks all tensors in SU(2)L representations. In particular, the doublet QL is broken in QL1 and QL2 (the names are chosen automatically, adding _1 _2 _3... to broken fields names). We rename then the broken left quarks to uL and dL:
toyModel.renameParticle("Q_L_1", "u_L");
toyModel.renameParticle("Q_L_2", "d_L");
Simulating a Higgs mechanism
Here we do not want to have a complete Higgs mechanism, juste to add mass terms ourselves. A Higgs mechanism can be implemented as well and one may found examples in the manual. This will in particular introduce a bit of CSL, how to create mathematical expressions.
First, we needed to simulate the rotation between B and W3 that gives the photon A and the Z boson as well as the linear combination of W1 and W2 giving W+:
Particle W_1 = toyModel.getParticle("W_1");
Particle W_2 = toyModel.getParticle("W_2");
Particle F_W_1 = W_1->getFieldStrength();
Particle F_W_2 = W_2->getFieldStrength();
Particle W = W_1->generateSimilar("W");
W->setSelfConjugate(false); // W is complex, not W_1
Index mu = MinkowskiIndex();
Index nu = MinkowskiIndex();
// W_1 goes to (W+ + W-) / sqrt(2)
toyModel.replace(W_1, (W(mu) + GetComplexConjugate(W(mu))) / sqrt_s(2));
// Same for the field strength F_W
toyModel.replace(F_W_1, (W({mu, nu}) + GetComplexConjugate(W({mu, nu}))) / sqrt_s(2));
// W_2 goes to i*(W+ - W-) / sqrt(2)
toyModel.replace(W_2, CSL_I * (W(mu) - GetComplexConjugate(W(mu))) / sqrt_s(2));
// Same for the field strength F_W
toyModel.replace(F_W_2, CSL_I * (W({mu, nu}) - GetComplexConjugate(W({mu, nu}))) / sqrt_s(2));
Expr thetaW = constant_s("thetaW");
Expr cW = cos_s(thetaW);
Expr sW = sin_s(thetaW);
// We give the rotation matrix explicitly here, between double curly braces
toyModel.rotateFields(
{"W_3", "B"},
{"Z" , "A"},
{{cW, sW},
{-sW, cW}}
);
Expr e = constant_s("e");
Expr g_Y = toyModel.getScalarCoupling("g_Y");
Expr g_L = toyModel.getScalarCoupling("g_L");
toyModel.replace(g_Y, e / cW);
toyModel.replace(g_L, e / sW);
There is here already a bit of work with expressions and tensors. It is on purpose, to help the reader to quickly get good habits on the symbolic manipulations in CSL. In particular, when we introduce W+, we have to get indices in the Minkowski space to use it with the Particle objects between parentheses, and curly braces if there is more than one index. We also just saw how to get symbolic functions, sums, products, the imaginary i, constants, and how to store expressions in variables.
Expr M_W = constant_s("M_W");
Expr m_u = constant_s("m_u");
Expr m_d = constant_s("m_d");
toyModel.addBosonicMass("W", M_W);
toyModel.addFermionicMass("u_L", "u_R", m_u);
toyModel.addFermionicMass("d_L", "d_R", m_d);
This is a good bit of work ! Now let's refresh the model to let MARTY do some simplifications in the Lagrangian and see what is inside our particle physics model:
toyModel.refresh();
cout << toyModel << endl;
Calculations
Feynman rules
This part is really simple. Feynman rules are calculated automatically when launching a computation if they have not already been calculated. One may still get them explicitly, displaying them on the terminal and showing the diagrams in GRAFED:
auto rules = toyModel.getFeynmanRules();
Display(rules); // Displays expressions in terminal
Show(rules); // Shows diagrams in the application
Amplitude
Let's now apply our model to a physical calculation, the uū→dd̅ transition. This process is interesting because all the vector bosons participate to the it. To calculate an amplitude, MARTY needs field insertions, and a perturbation theory order. Let's consider the tree-level to simplify the example. Keep in mind that the rest of the code would be the same for other insertions or at one-loop.
auto res = toyModel.computeAmplitude(
Order::TreeLevel,
{Incoming("u"), Incoming(AntiPart("u")),
Outgoing("d"), Outgoing(AntiPart("d"))}
);
We gave the order, the model, and insertions given by names, specifying each time
if the field is incoming or outgoing, and particle (default) or anti-particle.
The result is stored in res, that contains expressions (res.expressions)
and diagrams (res.diagrams). Expressions can be displayed in the
terminal with Display(res) and diagrams can be showed by
GRAFED with Show(res):
Display(res);
Show(res);
Squared amplitude
Then, let's create with MARTY a library containing the squared amplitude of this process as a numerical C++ function. First, we need to calculate this quantity:
Expr squared_ampl = toyModel.computeSquaredAmplitude(res);
The squared amplitude is a scalar symbolic expression. It is averaged over incoming spins.
Library generation
Now, let's place the library in a given location and compile it automatically with MARTY:
mty::Library myLib("toy", "<library-location-path>");
myLib.addFunction("squared_ampl", squared_ampl);
myLib.build();
One has to precise here mty:: because
it is an object that already exists in CSL (csl::Library), doing the same thing.
The specialized object in the namespace mty
simply adds some standard options for the generated library to compile properly.
It is recommended to always use the library in the mty
namespace.
If MARTY and/or CSL are not installed in a standard path, one must add the following lines before building the library:
myLib.addIPath("<marty-installation-path>/include");
myLib.addIPath("<csl-installation-path>/include");
myLib.addLPath("<marty-installation-path>/lib");
myLib.addLPath("<csl-installation-path>/lib");
Library usage
Now that the library is built, let's go in its folder and try to evaluate our
squared amplitude. By modifying the file script/example_toy.cpp
and typing make in the root directory of the library, one may directly
get a running program in bin/example_toy.x. For example:
#include "toy.h"
using namespace std;
using namespace toy;
int main() {
// Initial parameter definitions
double M_W = 80.379;
double e = sqrt(4 * M_PI * 1. / 137);
double m_u = 0;
double m_d = 0;
double thetaW = asin(sqrt(0.231));
// Kinematics left to the user
double E_CM = 1e2; // 100 GeV in center of mass
double theta = M_PI / 3; // Let's take an angle of pi / 3
double s = E_CM * E_CM;
double t = s / 2 * (1 - cos(theta));
double u = s / 2 * (1 + cos(theta));
double s_12, s_13, s_14, s_23, s_24, s_34;
s_12 = s_34 = s / 2;
s_13 = s_24 = -t / 2;
s_14 = s_23 = -u / 2;
param_t params;
params.M_W = M_W;
params.reg_prop = 0; // Regulator not needed here
params.e = e;
params.m_u = m_u;
params.m_d = m_d;
params.thetaW = thetaW;
params.s_12 = s_12;
params.s_13 = s_13;
params.s_14 = s_14;
params.s_23 = s_23;
params.s_24 = s_24;
params.s_34 = s_34;
// Final result
cout << "Squared amplitude: ";
cout << squared_ampl(params);
cout << endl;
return 0;
}
Finally, here is the result:
Squared amplitude = (0.201037,0)
One can see that the squared amplitude is real, and the final value will of course depend on the initial parameters. We built a model entirely from scratch, calculated an amplitude followed by its squared value implying even more simplifications, obtained its analytical value and generated automatically a numerical C++ library evaluating this analytical expression for a given set of initial parameters. All of this in less than a hundred lines of code.