What is xloops?

xloops was written to implement the methods developed in Mainz for calculating one- and two-loop Feynman diagrams. Whereas many other packages for loop calculations are often ad-hoc implementations computing numeric results for very specific processes or special cases (but quite efficiently so), the idea of xloops is to have a program that works for arbitrary massive diagrams to any tensor rank, and that can also return analytic results if possible. Eventually, xloops will be expanded into a package for the complete automatic calculation of processes in the Standard Model.

The roots of xloops go back to a Maple package for one-loop integrals released in 1993. Due to inherent problems with Maple as an implementation environment, the project was switched to C++ in 1999. One result of this was the creation of GiNaC, a C++ library for symbolic computation, that is now a project of its own.

Currently, xloops consists of a library of C++ and GiNaC functions for calculating arbitrary one-loop one-, two-, and three-point, and two-loop self-energy integrals of any tensor rank. All one-loop integrals and the divergent part of two-loop integrals are returned in symbolic form as an expansion in the dimensional regularization parameter ε. The finite part of two-loop integrals is transformed into a two-dimensional function that is numerically integrable.


  1. A one-loop calculation
    Here is a complete example program that calculates the divergent and finite parts of the general rank-1 one-loop two-point integral:
    #include <iostream>
    #include <ginac/ginac.h>
    #include <xloops/xloops.h>
    using std::cout; using std::endl;
    using namespace GiNaC;
    using namespace xloops;
    int main(void)
        // Declare symbols for momentum, masses, and parameters
        symbol q("q"), m1("m1"), m2("m2"), eps("eps"), rho("rho");
        // One-loop integral function
        ex C = OneLoop2Pt(1, 0, q, m1, m2, 1, 1, rho, eps);
        // Perform epsilon-expansion
        C = C.series(eps == 0, 1);
        // Extract the divergent and finite parts
        ex div = C.coeff(eps, -1).subs(rho == 0).normal();
        ex fin = C.coeff(eps, 0).subs(rho == 0);
        // Output the result
        cout << "Order eps^-1 is " << endl << div << endl;
        cout << "Order eps^0 is "  << endl << fin << endl;
        return 0;
  2. A two-loop calculation
    The following program first computes the divergent part of the scalar two-loop sunrise topology. Then it evaluates the finite part of the integral for a specific choice of the masses and the external momentum.
    #include <iostream>
    #include <ginac/ginac.h>
    #include <xloops/xloops.h>
    using std::cout; using std::endl;
    using namespace GiNaC;
    using namespace xloops;
    int main(void)
        // Declare symbols for momentum, masses, and parameters
        symbol q("q");
        symbol m1("m1"), m2("m2"), m3("m3");
        symbol rho("rho"), eps("eps");
        symbol x("x"), y("y");
        // Two-loop two-point integral function, sunrise topology
        ex C = TwoLoop2Pt4(0, 0, 0, 0, 0, q, lst(m1, m2, m3), lst(1, 1, 1), rho, eps, lst(x, y));
        // The result is a two-element list containing the analytic and numeric
        // parts; extract them
        ex num = C.op(0);
        ex ana = C.op(1);
        // Perform epsilon-expansion of the analytic part
        ana = ana.series(eps == 0, 1);
        // Extract the divergent (eps^-2 and eps^-1) and finite parts
        ex eps2 = collect_common_factors(ana.coeff(eps, -2).subs(rho == 0).normal());
        ex eps1 = collect_common_factors(ana.coeff(eps, -1).subs(rho == 0).normal());
        ex eps0 = ana.coeff(eps, 0);
        // Output the result for the divergent part
        cout << "eps^-2 = " << eps2 << endl;
        cout << "eps^-1 = " << eps1 << endl;
        // Substitute values for masses and momentum in the finite parts
        lst s;
        s = rho == 1E-15,
            m1 == 2,
            m2 == 4,
            m3 == 6,
            m4 == 8,
            m5 == 10,
            q == 20;
        eps0 = eps0.subs(s).evalf();
        fin = fin.subs(s);
        // Perform numerical integration using Vegas
        fin = VNumInt(fin, 10000, 10, 100000, 10, lst(x, y));
        // Output the result for the finite part, including the finite
        // analytic term
        cout << "finite = " << fin.op(0) + fin.op(2)*I + eps0 << endl;
        cout << "error real part = " << fin.op(1) << endl;
        cout << "error imag part = " << fin.op(3) << endl;
        return 0;