ColorFull -- a C++ library for calculations in SU(Nc) color space

ColorFull, a C++ package for treating QCD color structure, is presented. ColorFull, which utilizes the trace basis approach, is intended for interfacing with event generators, but can also be used as a stand-alone package for squaring QCD amplitudes, calculating interferences, and describing the effect of gluon emission and gluon exchange.


Introduction
The description of QCD color structure in the presence of many external colored partons is a field of increased importance. Some methods for performing automatic color summations of fully contracted vacuum bubbles, for example as implemented in FeynCalc [1], in the C program COLOR [2], or as presented in [3], have been around for a while, and recently a more flexible Mathematica package, ColorMath [4], allowing color structures with any number of open indices, has been published. Yet other general purpose event generator codes, such as MadGraph [5], have separate built in routines for dealing with the color structure.
In the present paper a stand-alone C++ code, ColorFull, designed for dealing with color contraction using color bases is presented 1 . ColorFull is written with interfacing to event generators in mind, and is currently interfaced to Herwig++ (2.7) [6,7], but can also be used as a stand-alone package for investigations in color space.
The intent of this paper is to convey the underlying idea of ColorFull. For full technical details we refer to the online Doxygen documentation 2 . To set the stage, a brief introduction to the QCD color space is given in section 2 and the trace basis approach is presented in section 3. Following this, some remarks about the computational strategy are made in section 4 and the key design features are presented in section 5, whereas examples of stand-alone usage are given in section 6, and the interface to Herwig++ is commented upon in section 7. The following section, section 8, describes the classes of ColorFull and code validation is discussed in section 9. Finally some concluding remarks are made in section 10.

Color space
Apart from four-gluon vertices for which the color structure can be rewritten in terms of (one-gluon contracted) triple-gluon vertices, the QCD Lagrangian contains quark-gluon vertices, where we follow the convention of reading the fully anti-symmetric structure constant indices in counter clock-wise order. The color structure of any amplitude, tree-level or beyond, pure QCD or not, can thus be expressed in terms of these objects alone. For observables we are -as QCD is confining -only interested in color summed/averaged quantities. Letting c 1 denote the color structure of the amplitude under consideration, we are thus interested in |c 1 | 2 where the scalar product 3 is given by summing over all external color indices, i.e., c 1 |c 2 = a 1 , a 2 , ... c * a 1 a 2 ... Clearly, in any QCD calculation, the color amplitudes, c 1 and c 2 , may be kept as they are, with color structure read off from the contributing Feynman diagrams. Alternatively -and this is likely to be the preferred solutions for more than a few partons -they may be decomposed into a color basis (spanning set), such as a trace basis [8][9][10][11][12][13][14][15][16], a color flow basis [19] or a multiplet basis [17,18].

Trace bases
One way of organizing calculations in color space is to use trace bases [8][9][10][11][12][13][14][15][16]. To see that this is always possible, we note that the triple-gluon vertex can be expressed as where T R is the normalization of the SU(N c ) generators, tr(t a t b ) = T R δ ab , commonly taken to be 1/2 or 1. Using this relation on every triple-gluon vertex in any amplitude results in general in a sum of products of (connected) traces over SU(N c ) generators and open quark-lines. More specifically, there is one open quark-line for every incoming quark/outgoing anti-quark and outgoing quark/incoming anti-quark. (Note that, from a color perspective outgoing antiquarks are equivalent to incoming quarks; we will here refer to them collectively as quarks. Similarly outgoing quarks are equivalent to incoming anti-quarks, and will be referred to as anti-quarks.) To further simplify the color structure, we may contract every internal gluon propagator (which now connects quark-lines) using the Fierz (completeness) relation From this we see that every amplitude in QCD, at any order, may be expressed as a sum of products of open and closed quark-lines, where all gluon indices are external indices. The set of all such products of open quark-lines can thus be used as a spanning set for any QCD process. We here refer to this type of basis as a trace basis, although we remark that when the number of gluons, N g , plus the number of qq-pairs, N q , exceeds N c , this spanning set is overcomplete, and hence strictly speaking not a basis. For N g + N q ≤ N c , the bases are not overcomplete, but they are still non-orthogonal, having non-diagonal scalar products being suppressed by powers of N c . Only in the N c → ∞ limit, are these bases minimal and orthogonal.
To all orders in perturbation theory in the N c → ∞ limit, one can prove that the number of basis vectors can be found using the recursion relation [17] N vec [N q , , (3.4) where N vec [N q , 0] = N q ! , N vec [N q , 1] = N q N q ! . (3.5) In the gluon-only case, at tree-level, the only color structures that can appear are traces over generators, meaning that a general tree-level gluon amplitude can be decomposed as M(g 1 , g 2 , . . . , g n ) = σ∈S Ng −1 tr(t g 1 t gσ 2 . . . t gσ n )A(σ). (3.6) That only fully connected color structures enter in tree-level gluon amplitudes can easily be understood from the decomposition of Feynman diagrams into basis vectors; upon application of eq. (3.1) all external gluons remain attached to a quark-line, andwhile contracting internal gluons using the Fierz identity, eq. (3.2) -this remains true, as the color suppressed terms cancel each other out. The same cancellation appears for gluon exchange between a quark and a gluon, meaning that also tree-level color structures for one qq-pair and N g gluons must be of the "fully connected" form of a trace that has been cut open, an open quark-line, M(q 1 , g 3 , . . . , g n , q 2 ) = σ∈S Ng (t gσ 1 t gσ 2 . . . t gσ n ) q 1 q 2 A(σ), (3.7) i.e., only the first two basis vectors in eq. (3.3) are needed. However, when the Fierz identity is applied directly to a gluon exchange between quarks, as in eq. (3.2), both terms do appear, and color structures with up to N q disconnected quark-lines may appear even at tree-level. Starting from a tree-level color structure and exchanging a gluon between two partons may split off a disconnected color structure. Thus, counting to l g additional gluon exchanges (on top of a tree-level diagram), the color structures can not consist of more than max(1, N q ) + l g open and closed quark-lines. In general, when any Feynman diagram is decomposed into the trace basis, there can be at most N q + ⌊N g /2⌋ quark-lines, since all gluons may be disconnected from the quarks, but no gluon can stand alone in a trace.
For NLO color structures having a quark-loop in the Feynman diagram, the quark-loop is necessarily connected to the remaining color structure via a gluon exchange, unless the Feynman diagram consists of only gluons attached to one quark-loop, in which case the color structure is trivial. For each color structure connected to the quark-loop, contracting the connecting gluon may give rise to a disconnected color structure -if the connecting gluon goes to a quark.
There are at most N q gluons of this type, giving, rise to at most N q + 1 quark-lines in the fully decomposed color structure, i.e., for processes with quarks, one more than at treelevel. (For purely gluonic processes, no additional trace will appear.) The same argument applies if the Feynman diagram contains more than one closed quark-loop, giving (also in this case) up to max(N q , 1) + l q color structures for l q quark-loops.
Finally, we note that, if, on top of the quark-loops, there are also l g loops from gluon exchange, these may contribute with up to one additional quark-line each, completing the argument that for n l = l q + l g loops, a general color structure in pure QCD can be written in terms of at most min[max(1, N q ) + n l , N q + ⌊N g /2⌋] closed and open quark-lines.

Computational strategy
The ColorFull strategy for treating color space is based on the above observations, i.e., • For given external particles (quarks, anti-quarks and gluons), we may always decompose the color space into a linear combination of closed and open quark-lines, as described in section 3. (Other color bases may be expressed in terms of linear combinations of these terms.) • We can always evaluate scalar products by first replacing triple-gluon vertices using eq. (3.1), then removing all gluon propagators using eq. (3.2), and finally contracting the remaining product of qq-delta functions.
The above outlined procedure for calculating scalar products has the advantages of being conceptually easy, and of covering all contractions occurring in QCD. It has the disadvantage of potentially giving a very large number of terms, since every replacement of a structure using eq. (3.1) may double the number of terms, and similarly, so does naive application of the Fierz identity, eq. (3.2).
In order to mitigate this potential explosion of terms upon index contraction, the ColorFull procedure for scalar product contraction is: 1. Contract all quark-ends, giving a sum of products of closed quark-lines.
2. On the individual quark-lines, contract all neighboring gluons (giving a factor C F ) and all next to neighboring gluons, giving a factor −T R /N c each, as only the color suppressed term in the Fierz identity, eq. (3.2), survives. Also, contract traces of two gluons, tr(t g 1 t g 2 ) = δ g 1 g 2 /T R .
3. Look for gluons to contract within the quark-lines. Each such contraction may give rise to two terms, but at least the new traces tend to be shorter.

Look for gluon index contractions between the quark-lines.
In this way all color indices can be iteratively contracted. To further speed up calculations, ColorFull can use memoization to save calculated topologies, meaning that contractions that differ only by a relabeling of indices are performed only once. This significantly speeds up the calculations.

ColorFull at a glance
ColorFull is designed to handle the contraction of QCD color indices, to decompose the QCD color space using (trace) basis vectors, and to describe the effect of gluon exchange and gluon emission. In particular, ColorFull can, for arbitrary N c : • Square any QCD amplitude and calculate any interference term.
• Create a trace basis for any number of quarks and gluons, and to arbitrary order.
• Read in and write out color bases, including non-trace bases.
• Calculate scalar product matrices, i.e., the matrices of scalar products between the basis vectors, write these out and read them in again.
• Describe the effect of gluon exchange, including calculating the color soft anomalous dimension matrices.
• Describe the effect of gluon emission.
ColorFull can also be interfaced to Herwig++ (≥ 2.7) via Matchbox [20] and can thus easily be used for event generation with Herwig++.

Class Functionality
Monomial Class for containing a single term of form Nc pow Nc × TR pow TR × CF pow CF × int part × cnum part and associated functions.

Polynomial
Class for containing a sum of Monomials and associated functions.

Poly vec
Class for containing a vector of Polynomials and associated functions.

Poly matr
Class for containing a matrix of Polynomials and associated functions.

Quark line
Class for containing a single closed or open quark-line multiplying a Polynomial and associated functions. For example a term of the form tr[t g 1 t g 2 t g 3 t g 4 ] ∼ (1,2,3,4) or of the form (t g 3 t g 4 ) q 1 q 2 ∼ {1,3,4,2} times a Polynomial.

Col str
Class for containing a product of Quark lines multiplying a Polynomial, and associated functions. For example a polynomial (N 2 c − 1) may multiply the Quark lines, δ q 1 q 2 δ g 1 g 4 , in total giving the Col str

Col functions
Library class containing functions for performing numerical evaluations, taking the leading N c limit, evaluating scalar products, and describing gluon emission and exchange.

Col basis
Base class for all basis classes, see below.
Trace type basis Base class for Trace basis and Tree level gluon basis.

Trace basis
Class for creating and using trace bases. Perhaps the most important basis class.

Tree level gluon basis
Class for bases needed for tree-level gluon calculations, where each basis vector is a single trace plus an implicit complex conjugate.

Orthogonal basis
Class for utilizing the benefits of orthogonal bases such as multiplet bases [17,18].

Technical overview
This section is intended to give an overview of ColorFull. For examples, the reader is referred to section 6, for overviews of the classes to section 8.2-8.3, and for details, the tables in the appendix A, as well as the online Doxygen documentation. As contraction of indices gives rise to polynomials in N c , T R and C F , ColorFull by necessity needs minimal classes for dealing with such polynomials. This is implemented in the classes Monomial (a single term in a polynomial 5 ), Polynomial (a sum of Monomials), Poly vec (a vector of Polynomials) and Poly matr (a matrix of Polynomials).
For the color structure itself, ColorFull uses the class Quark line for treating an individual closed or open quark-line, a class Col str for treating a product of Quark lines and a class Col amp for treating a general color amplitudes, i.e., a linear combination of Col strs.
For performing numerical evaluations, taking the leading N c limit, evaluating scalar products, and describing gluon emission and exchange, ColorFull has a library class Col functions.
Finally, ColorFull has classes for describing color bases. The classes intended for the user, Trace basis, Tree level gluon basis and Orthogonal basis, are derived from the base class Col basis(in the case of the first two via the class Trace type basis).
All ColorFull classes are listed in table 1 according to dependencies, meaning that each class only depends on classes listed above it. The next section will give an introduction to using ColorFull in stand-alone mode.

Stand-alone usage
ColorFull is mainly designed to deal with relatively large color spaces where it is advantageous to use bases in which color structures coming from Feynman diagrams -or alternative recursive strategies -are decomposed.
For trace bases, the bases may be automatically created by ColorFull. For example, a basis for 1 qq-pair, 3 gluons and 0 loops (in pure QCD) can be created using Trace basis MyBasis(1,3,0);. (6.1) The last argument is provided to avoid carrying around basis vectors for which the kinematic factors vanish to a certain order in perturbation theory. It can be skipped upon which an all order basis vector is created.
To view the resulting basis, it can be written out to a stream (cout) or to a file, either with default filename (no argument) or a user supplied name, cout << MyBasis; MyBasis.write out Col basis("ColorResults/MyBasis"); (6.2) 5 Throughout the ColorFull documentation terms of form constant×N a c T b R C c F will be referred to as monomials, despite the possible occurrence of negative powers. Similarly sums of such terms will be referred to as polynomials. where, for example, If no argument is supplied to MyBasis.write out Col basis(), the basis is written out to a file with a default filename, in this case CF TB q 1 g 3. The bases can thus be written out and saved for future purposes. For reading in a basis MyBasis.read in Col basis("path/to/filename") can be used. The option of reading in bases is particularly useful for the Orthogonal basis class. Presently ColorFull can not automatically create orthogonal multiplet bases, such as in [17], but externally created (orthogonal) bases, with basis vectors expressed in terms of sums of products of traces, may be read in and used. For the Orthogonal basis class, the orthogonality is then utilized to significantly speed up the the calculation of scalar products. For a non-orthogonal basis it is necessary to evaluate all scalar products between all basis vectors, which can be done as MyBasis.scalar product matrix();. (6.5) Polynomial and double versions of the scalar product matrix are then calculated and saved in the member variables P spm and d spm. For larger bases, the evaluation can be sped up by only calculating the numerical version using MyBasis.scalar product matrix num(). The content of P spm and d spm can be explored by using the << operator, but it may also be saved to a file using the Col basis member functions write out P spm() and write out d spm() which write out the result in ColorFull and Mathematica, readable format, by default to a file with a default filename in the directory ColorResults. Alternatively the user may supply a filename as argument.
In several contexts, such as parton showers, cancellation of infrared singularities in NLO calculations, and recursive methods for calculating amplitudes, it is of interest to know the effect of gluon emission on a color structure. This can be calculated by using the Col functions member function emit gluon Col fun.emit gluon(Ca1,3,6); The sign conventions here and elsewhere are such that every gluon inserted after the emitter in the quark-line comes with a plus sign and every gluon inserted before comes with a minus sign. In a basis-decomposed calculation, one would be interested in this result decomposed in the large basis required for one additional gluon. Having a trace basis for this larger vector space, this decomposition can be calculated using Trace basis LargerBasis(nq,ng+1); (6.10) LargerBasis.new vector numbers(Cs, emitter); (6.11) where Cs is the Col str (basis vector in the smaller trace basis) from which the parton emitter emits a new gluon. Similarly the effect of gluon exchange on a color structure is of interest. ColorFull offers several functions for dealing with this. Starting with an amplitude, the new amplitude after gluon exchange can be obtained as In some situations, such as soft gluon resummation, it is also useful to calculate the soft anomalous dimension matrix, i.e., to have the result of gluon exchange on any basis vector decomposed into the basis. This can be computed automatically using the Col basis member function color gamma. The result is contained in a Poly matr. For this we need the full basis since color structures not present at LO will appear at NLO, etc. Coding In this way all the soft anomalous dimension matrices needed in [21] can easily be recalculated. While the number of colors in QCD is three, ColorFull can deal with any N c , both algebraically and numerically. Numerical evaluation is handled by the Col functions class, using the (private) member variables Nc, TR and CF. Thus also T R and C F can be changed independently. The reason for keeping C F as a parameter technically independent of N c is that this allows for keeping the color suppressed part of C F , −T R /N c , in a consistent way. This choice has proved useful for accounting for sub-leading N c effects in several phenomenological studies [22,23].
To numerically evaluate a Monomial, Polynomial, Poly vec or Poly matr, the Col functions member functions double num are used, for example we may want the N c = 3 version of a scalar product Col fun.double num(Col fun.scalar product(Ca1,Ca1)); (6.16) For comparison, it is of interest to evaluate (squared) amplitudes in the limit N c → ∞. For taking the leading N c limit of any of the polynomial classes Polynomial, Poly vec or Poly matr, the Col functions member function(s) leading can be used. For example we may take the leading N c limit of the matrix in eq. (6.15), For an expression containing C F , the color suppressed term −T R /N c can be kept in numerical evaluation by setting full CF to true, Col fun.set full CF(true).

Interfacing to Herwig via Matchbox
ColorFull can also be interfaced to Herwig++ (≥2.7) [7] via the Matchbox component [20], and can be used to treat the hard interaction, as in for example [24], as well as the parton shower itself [22]. When linked to Herwig++, ColorFull is hooked into the boost linear algebra package, enabling a very efficient treatment of numerical linear algebra. From the next major release of Herwig++, ColorFull will be directly shipped with the Herwig++, code.

Operators
In order to simplify calculations and increase readability, ColorFull defines a standard set of operators, listed in table 2. This includes the arithmetic operators +, -and * as well as the comparison operators, == and !=, and the stream operator << . The comparison operators, == and !=, work by comparing data entries term by term, both for the polynomial classes and for the color carrying classes. Thus, for example, two Polynomials only differing by the order of terms are not considered equal.

Classes for Polynomials
The result of color index contraction can always be written as a sum of terms of the form N c a T R b C F c × constant where we allow for negative integers a, b, c. For the purpose of treating contracted color structures, ColorFull has a minimalistic set of classes for basic manipulation of color factors arising when contracting color indices. One such term is defined as a Monomial, and a sum of such terms as a Polynomial. To decompose vectors in a color space we also need a vector of Polynomials, contained in a Poly vec, and for a scalar product matrix (or soft anomalous dimension matrix) we need a matrix of polynomials, a Poly matr. Apart from the Monomial class, these classes have the actual polynomial information contained in uncapitalized typedefs carrying the same name as the class in question. There are thus polynomial, poly vec, and poly matr typedefs. For manipulating the polynomial classes, the operators listed in table 2 may be used. For example, we may multiply a Polynomial and a Monomial.

Monomial
The simplest class for containing contracted color information is a Monomial The exponents pow Nc, pow TR and pow CF, as well as int part are of type int and cnum part is a cnum, an std::complex<double>. These member variables carry the information about the actual monomial. By default, using the standard constructor, Monomial(), a Monomial is set to 1, by letting cnum part=1, int part=1 and setting all powers to 0 6 . A Monomial with different integer part can be obtained by using the constructor taking an int or as argument, and a general Monomial can be constructed using the string constructor (see table 3).

Polynomial
A sum of Monomials is contained in the class Polynomial where, technically, the polynomial information is stored in the public member poly, of type polynomial=std::vector<Monomial>. ColorFull is not designed for manipulating polynomials but Polynomial nevertheless has a minimalistic simplify() function which collects terms with equal power of N c , T R and C F . It also has a remove CF() function for replacing C F with T R N c − T R /N c . Details about the Polynomial class can be found in table 4 in appendix A.

Poly vec
For dealing with vectors of Polynomials, ColorFull has a class Poly vec Poly vec is a container class for functions acting on vectors of Polynomials whereas the actual information is stored in the poly vec=typedef std::vector<Polynomial> data member, pv. The member functions of Poly vec, which include simplify() and remove CF(), can be found in table 5 in appendix A. 6 A default Monomial is thus the neutral element under multiplication, rather than under addition.

Poly matr
Finally, there is a class for treating matrices of Polynomials, a Poly matr, which stores the matrix information as a vector of vectors of Polynomials, a poly matr=typedef std::vector<Poly vec>. Again, details can be found in the appendix, in table 6.

Classes containing the color structure
This section describes the building blocks used by ColorFull to treat the color structure. All classes used to carry the color structure, Quark line, Col str, Col amp, including the basis classes Col basis -from which Trace type basis, Trace basis, Tree level gluon basis and Orthogonal basis are derived -have the property that the actual information of the color structure is contained in a type with corresponding name, whereas the class acts as a container for related functions. Thus, for example, the color amplitude information in a Col amp is contained in a col amp variable.

Quark line
The most basic color carrying class is a Quark line. Loosely speaking, a Quark line is a quark line times a Polynomial, or the open Quark line Quark lines may be created using the Quark line(std::string) constructor. Closed quark-lines are denoted by standard parenthesis, whereas curly brackets represent open quark-lines. For example, we may write

Col str
A general color amplitude consists not only of one quark-line, but of a linear combination of products of Quark lines. One term in this linear combination, a Polynomial times a product of closed and open quark-lines is contained in a Col str, Here the actual information about the color structure is carried by a col str, (technically a vector of Quark lines), for example we may thus have Like the other color classes, Col strs may be created using a string constructor where we note that the col str is written inside square brackets. The Polynomial should multiply the whole col str, rather than individual quark lines.

Col amp
The linear combination of Col strs which in general is needed to keep track of a color structure is contained in a Col amp. In order to contain scalar results, i.e., terms not containing any color structure (which arise in the process of color contraction) a Col amp also contains a Polynomial part, Scalar, .
For example, we may have the Col amp where Scalar is 0, and the Polynomials multiplying the Col strs (1,2)(3,4) and

Col functions -a library class
Col functions is a library class containing functions which cannot, in a natural way, be attributed to one class, or functions which act on many classes and are therefore conveniently collected into one class. In particular, Col functions contains classes for evaluating scalar products, for numerical evaluation, for taking the leading N c limit, and for describing the effect of gluon emission or gluon exchange.

Numerical evaluation
Col functions is the class which carries numerical values of N c , T R and C F , stored in the private members Nc, TR and CF respectively. To numerically evaluate a Monomial to a double or a complex number (cnum) the function cnum num or double num is used The syntax for numerical evaluation of Polynomials, Poly vecs and Poly matrs is identical.

Leading N c evaluation
Col functions also contains functions for taking the leading N c limit of the classes Polynomial, Poly vec and Poly matr. The leading N c terms can be evaluated in two different ways, either by taking the strict N c → ∞ limit and dropping all color suppressed terms (default), or by keeping the color suppressed terms in C F = T R (N 2 c − 1)/N c , while dropping other color suppressed terms. For keeping the full C F in numerical evaluations the member variable full CF must be set to true using set full CF(true). Col fun.scalar product(Cs1,Cs2);

Scalar products
in both cases resulting in T R N c C 3 F .

Gluon exchange and gluon emission
For soft gluon resummation, for NLO calculations, and for the cancellation of real emissions and virtual corrections in the soft limit, the color structure associated with gluon exchange plays an important role. The Col functions class therefore has functions for describing the effect of gluon exchange on a Col str, or on a Col amp. For example, we may exchange a gluon between parton 3 and 6 in Cs1 Col fun.exchange gluon(Cs1,3,6); To understand the signs we note that each time a gluon is inserted before the emitter on a (closed or open) quark-line there is a minus sign, and each time a gluon is inserted after, there is a plus sign. We also note that in this case, the result of gluon exchange on a single color structure gave rise to a linear combination of four color structures, the maximal possible number of color structures from a single col str [16].
In the context of gluon exchange we also remark that Col functions can calculate the "color correlator", c|T i · T j |c arising when coherently emitting a gluon from parton i and j in an amplitude |c , or when exchanging a gluon between the partons i and j in |c . For example we can calculate the color correlator for exchanging a gluon between parton 1 and 4 in Ca1, Col fun.color correlator(Ca1,1,4);

Classes for color bases
Although ColorFull may perform calculations with individual Quark lines, Col strs and Col amps, the intended usage is via the color basis classes Trace basis (in particular), Tree level gluon basis and Orthogonal basis.
As these classes share much of the most important functionality, they all inherit from one base class, Col basis. The Trace basis and Tree level gluon basis classes inherit from Col basis via Trace type basis, whereas Orthogonal basis inherits directly from Col basis.

Col basis
The Col basis class has a col basis member variable cb for containing the basis vectors, The most important Col basis member function is probably the scalar product matrix() function which calculates the matrix of scalar products between the basis vectors. In its default from, this function uses memoization, as this speeds up the calculations, but this may be circumvented in using the scalar product matrix no mem() version.
In order not to have to calculate scalar product matrices over and over again Col basis also has functions for reading in and writing out scalar product matrices, both in numerical and polynomial form.
Col basis also contains functions for reading in and writing out the basis itself. This is essential for the Orthogonal basis class which cannot (presently) construct the orthogonal bases, but may read them in.
The decomposition of color amplitudes into bases is done with the (virtual) decompose(Col amp) member function, which does the decomposition by exploring the coefficients in front of traces and products of traces for Trace basis and Tree level gluon basis and by evaluating scalar products for the Orthogonal basis case.
Another important function is the color gamma function which (using decompose) calculates the soft anomalous dimension matrix, i.e., calculates the matrix describing the effect of gluon exchange between two partons. The result is returned as a Poly matr where component i, j gives the amplitude for ending up in vector i, if starting in vector j, see also section section 6. A list of public members and functions for Col basis is found in table 11.

Trace type basis
Trace type basis is a small helper class for keeping track of functions which are similar for Trace basis and Tree level gluon basis. It inherits from Col basis and is inherited from by Trace basis and Tree level gluon basis. Most importantly, this is where decompose is implemented for these two classes, see table 12.

Trace basis
Although the observation that each color amplitude may be decomposed into products of open and closed quark-lines is a guiding principle for ColorFull, the Trace basis class itself is rather small, containing mainly functions for creating bases.
A trace basis is created by first contracting all qq-pairs in all N q ! ways, and then attaching gluons to these closed quark-lines, and to additional closed quark-lines, such that at least two gluons are attached to each closed quark-line.
This can be done either directly using a constructor At tree-level, in pure QCD, the last two basis vectors vanish. To create a basis which is only valid up to order n loop in pure QCD, we may use the create basis(n quark,n gluon, n loop) member function.

Tree level gluon basis
In the case of gluon-only color structures, charge conjugation implies that each trace must appear with its conjugate, in linear combinations of the form tr[t g1 t g2 ...t gn ] + (−1) Ng tr[t gn ...t g2 t g1 ]. In the trace type basis class Tree level gluon basis this is used to reduce the number of basis vectors and speed up calculations.

Orthogonal basis
ColorFull can not -in its present form -automatically create multiplet bases. However, orthogonal bases may be read in using the (Col basis) function read in Col basis(std::string) member function. For dealing with orthogonal bases, ColorFull offers special functions for calculating the matrix of scalar products and decomposing vectors. The syntax for basis reading is the same as for basis writing. Bases thus appear much as as in eq. (8.23), see table 11.

Validation
For a code with order 10 000 lines, validation is essential. For this reason ColorFull is continuously validated using a test suite, which aims at testing all the various components. The applied tests starts with checking basic functions for reading in and writing out files, and dealing with polynomials. After this, the creation of bases is tested, and scalar products are tested by changing the order of index contraction, and by switching on and off memoization. The scalar product matrices have further been tested against ColorMath [4] and, for processes occurring in the context of [24], also against another private Mathematica code. The functions describing gluon emission and gluon exchange are cross-checked against each other.

Conclusion and outlook
ColorFull, a C++ stand-alone QCD color algebra package, designed for interfacing with event generators, has been presented.
ColorFull is based on trace bases, which can automatically be created, and color contraction is performed by repeated usage of the Fierz identity. Employing these bases, one can in principle describe any QCD process. In reality, the scalar product matrices, which may be calculated once and for all, become hard to manage for more than approximately 8 gluons plus qq-pairs. This limitation is inherent for the trace bases, since they are nonorthogonal, and for this reason ColorFull is written to be able to load and use orthogonal (multiplet) bases.
ColorFull does, however, not -in its present form -perform index contraction in terms of group invariants as described in for example [18,25]. Extending ColorFull to inherently construct orthogonal multiplet bases and efficiently perform index contraction using 3j and 6j coefficient may speed up the treatment of QCD color space very significantly.

Acknowledgments
Simon Plätzer is thanked for interfacing ColorFull to Herwig++ and for writing the build system. I also want to thank Simon Plätzer and Johan Grönqvist for many valuable discussions on the organization of this code.  Constructor using a double. The cnum part member is set to contain the value.

Monomial(int num)
Constructor using an int. The int part member is set to contain the value.
Monomial(std::string str) Constructor taking a string as argument. The argument should be of the form in for example "-(20*TR^5)/Nc" or "-20 TR^(5)/Nc" or "20 / TR^(-5)Nc^(1) CF^(3)". Note: All spaces and all * are ignored, except in *(-1) and *-1, which are understood as *(-1). Everything standing after / is divided out, whereas everything standing before / is multiplied with. Parentheses are ignored unless they appear in powers, directly afterˆ. No spaces are allowed inside the powers. If the string contains no information or is empty, the Monomial is set to 1, i.e., pow TR = pow Nc = pow CF = 0, int part = 1, cnum part = 1.0.    An empty polynomial is defined as 1, to get 0, multiply with 0.

Polynomial(double dnum)
Constructor allowing setting the Polynomial using a double. The Polynomial gets one Monomial where the real part of cnum part equals dnum.

Polynomial(int num)
Constructor allowing setting the Polynomial using an int.

Data member Content
poly vec pv where poly vec is: typedef std::vector<Polynomial> To actually contain the vector of Polynomials.

Poly vec(poly vec poly v)
Makes a Poly vec of a poly vec.

Member functions Effect void append(Polynomial Poly)
Appends a Polynomial to data member pv.
const Polynomial& at(int i) const Returns the Polynomial at place i.

Polynomial& at(int i)
Returns the Polynomial at place i. Writes out the vector to the file filename. Table 6. Public members and functions of the class Poly matr.

Data member Content
poly matr pm where poly matr is: typedef std::vector<Poly vec> To actually contain the matrix of Polynomials.

Constructor Effect
Poly matr() Default constructor, leaves pm empty.

Member functions Effect
Poly vec& at(int i) Returns the Poly vec at place i.  To actually contain the color information, in order {q, g 1 , g 2 , ...g n , q} or (g 1 , g 2 , ...g n ).

bool open
Is the string open, with a quark in the beginning and an antiquark in the end, or not?

Polynomial Poly
Polynomial factor, multiplying the quark line.

Constructors Effect
Quark line( ) Default constructor, leaves ql empty.

Member functions Effect
Quark line after(int j) const Returns a Quark line where the ql member is changed to contain only partons after place j.
void append(int p) Appends parton p to the Quark line.
void append(std::vector<int> in ql) Appends a whole quark line to the Quark line.

int at(int j) const
Returns the parton at place j.
For closed quark lines j may be between -size and 2*size.

Quark line before(int j) const
Returns a Quark line where the ql member is changed to contain only partons before place j.  Erases the parton located at place. std::string find kind(int p) const Finds out if a parton is a quark, anti-quark or gluon, returns "q", "qbar" or "g" respectively. This function does not loop over all partons, but assumes that the parton is a gluon if it is in a closed Function for writing out the Col str to a file with name filename. Table 9. Public members and functions of the class Col amp.

Data members Content
col amp ca where col amp is: typedef std::vector<Col str> To actually contain the information about the Col strs, ca=Cs1+Cs2+Cs3+... .

Polynomial Scalar
Scalar is a Polynomial for collecting color factors appearing when the color structure has been fully contracted. The full color amplitude is Scalar+Cs1+Cs2+Cs3.... Scalar should thus be non-zero only if all indices can be contracted.

Col amp(Col str Cs)
Constructor converting a Col str to a Col amp.
Col amp(const std::string str) Constructor taking a string as argument.
The string should be of the form Polynomial1*col str1+Polynomial2*col str2,  Function for writing out the Col amp to a file with name filename. Table 10.
Some (private) data members and the public functions of the library class Col functions.

double CF
The value of C F = T R (N 2 c − 1)/N c , changed by the set CF function. Note that CF can be changed independently of Nc. double Nc The number of colors, used in numerical results, changed by the set Nc function.

double TR
The trace convention tr(t a t a )=TR (no sum), the normalization of the SU(N c ) generators, to be used in numerical evaluation. This value can be changed by the set TR function.

Constructor Effect
Col functions() Default constructor, sets Nc=3, TR=0.5, CF=4.0/3.0 and full CF=false. Polynomial color correlator(const Col amp Ca, int i, int j) const Calculates c|T i T j |c , the "color correlator" relevant for coherent gluon emission from parton i and parton j, or gluon exchange between i and j. The Ca thus corresponds to |c , and i and j are the partons involved in the emission (exchange). Function for emitting a gluon from a Col str. When the gluon is inserted before the emitter in a Quark line, the amplitude comes with a minus sign.

Col amp emit gluon(const Col amp & Ca in, int emitter, int g new) const
Function for emitting a gluon from a Col amp. When the gluon is inserted before the emitter in a Quark line, the amplitude comes with a minus sign.
Col amp exchange gluon(const Col str& Cs, int p1, int p2) const Function for exchanging a gluon between the partons p1 and p2 in the Col str Cs. When the gluon is inserted before the emitter in a Quark line, the amplitude comes with a minus sign.
Col amp exchange gluon(const Col amp & Ca, int p1, int p2) const Function for exchanging a gluon between two partons p1 and p2 in the Col amp Ca. When the gluon is inserted before the emitter in a Quark line, the amplitude comes with a minus sign.

double get CF() const
Returns the value of CF.
double get Nc() const Returns the number of colors.

double get TR() const
Returns the normalization of the generators, tr(t a t b ) = TRδ ab .

bool get full CF() const
Returns true if full CF is true and false otherwise.

int factorial(int i) const
The factorial of an int, 0! is defined as 1.  Poly matr leading P spm To contain the Polynomial version of the leading part of the scalar product matrix.

int ng
The number of gluons, initially set to 0 and (possibly) changed with create basis, or while reading in the basis.

int nq
The number of qq-pairs, initially set to 0 and (possibly) changed with create basis, or while reading in the basis.

Poly matr P spm
To contain the Polynomial version of the scalar product matrix.

Constructor and destructor Effect
Col basis() Default constructor, sets nq=ng=0, and the private data members trace basis = tree level gluon basis = orthogonal basis = false. virtual ∼Col basis() Destructor.

Member functions Effect
Col amp& at(const int i) Returns the Col amp (i.e. basis vector) at place i. Each type of color basis has to implement a function for decomposing an amplitude in the color basis.

bool empty() const
Is the col basis empty?
Col amp exchange gluon( uint vec, int p1, int p2) Function for finding the resulting Col amp after exchanging a gluon between parton p1 and parton p2 in the basis vector vec. bool is Orthogonal basis() const Is it an Orthogonal basis?
bool is Trace basis() Is it a Trace basis?
bool is Tree level gluon basis() Is it a Tree level gluon basis?
void leading scalar product matrix() Finds the leading Nc scalar product matrices, leading P spm and leading d spm. If the polynomial scalar product matrix, P spm has been calculated, P spm is used, otherwise P spm is first calculated and the leading Nc limit is then taken of P spm. int n gluon check() const Returns the number of gluons in the Col basis after checking that each Col str in each Col amp has the same number of gluons. Function for calculating scalar products algebraically using the basis and the scalar product matrix (Poly matr P spm) in the basis.
(Does add implicit conjugated part for Tree level gluon basis, as these terms are contained in the matrix of scalar products.) void scalar product matrix() Function for calculating the scalar products matrix. This function loops over all basis vectors and stores the value of the scalar product between basis vector i and basis vector j in the i,j -entry in P spm and d spm. The calculation is done using memoization. The symmetry of the scalar product matrix is not used for the calculation, instead it is checked that the resulting matrix is indeed symmetric.
void scalar product matrix num() Returns a standard filename, used for writing out scalar product matrices. If leading is true, " l" is appended to the filename. If "poly" is true "P " is added to the filename, and if it is false "d ", as in double, is added to the filename.

Data member (protected) Content int max ql
The maximal number of quark-lines allowed in the basis. This is used for constructing bases that only are valid up to a certain order in QCD, such that unused information need not be carried around. This function is intended for tree-level processes with at most 2 qq-pairs. It finds the new vector numbers in the basis for N q +N g +1 partons after radiating a new gluon from the parton emitter. This function does not actually use the cb, but only calculates the basis vector number, which makes it much quicker than the general version. The old vector has number old num, and there were, before emission, N q quarks (+ N q antiquarks) and N g − 1 gluons. For emission from a quark or an antiquark there is only one resulting color structure, and -1 is returned in the place of the absent color structure. The second vector, where the new gluon is inserted before the emitter comes with a minus sign in the new total amplitude. Table 13. Public member functions implemented in the class Trace basis. As Trace basis inherits from Trace type basis all the public members of Trace type basis and Col basis are also available, see table 12 and table 11.

Constructors Effect
Trace basis() Default constructor, calls the Trace type basis constructor and sets nq = ng = 0.
Trace basis(int n quark, int n gluon) Constructor for creating a trace basis for n quark qq-pairs and n gluon gluons.
Trace basis(int n quark, int n gluon, int n loop) Constructor for creating a trace basis for n quark qq-pairs and n gluon gluons, keeping only those color structures that can appear to order n loop in pure QCD. (Note: For electroweak interactions more color structures may be needed.)

Member functions Effect
void create basis(int n q, int n g) Creates a trace basis with basis vectors saved in the cb member. Keeps all possible basis vectors, i.e., the basis is valid to any order in perturbation theory.
void create basis(int n q, int n g, int n loop) Creates a trace basis with basis vectors saved in the cb member. Keeps only basis vectors consisting of at most n q plus n loop quark lines. Table 14.
Public member functions of the class Tree level gluon basis. As Tree level gluon basis inherits from Trace type basis which inherits from Col basis, all the public members of Col basis and Trace type basis are also available, see table 11 and table 12.

Constructors Effect
Tree level gluon basis() Default constructor.
Tree level gluon basis(int n g) Constructor for creating a tree-level gluon basis with ng gluons.

Data members Content dvec diagonal d spm
To contain information about scalar products as a dvec, i.e., entry i is the square of vector i.

Poly vec diagonal P spm
To contain information about scalar products as a Poly vec, i.e., entry i is the square of vector i.

Constructor Effect
Orthogonal basis():Col basis() Default constructor, puts private variable orthogonal basis =true and calls the constructor of Col basis.

Member functions Effect
Poly vec decompose(const Col amp & Ca) The decomposition of a Col amp in an orthogonal basis is done by calculating scalar products and dividing out the norm. The norm is evaluated numerically.
void diagonal scalar product matrix( bool save P diagonal spm, bool save d diagonal spm, bool use mem) Calculates the diagonal entries in the scalar product matrix, and (depending on arguments), saves them to the member variables diagonal P spm and diagonal d spm. This function is used by the Orthogonal basis version of scalar product matrix. std::string diagonal spm file name( const bool leading, const bool poly) const Creates a default filename for writing out diagonal scalar products. The boolean variable leading should be true if the name is for a leading Nc variable. The filename is then modified accordingly.

Polynomial scalar product(const Col amp & Ca1, const Col amp & Ca2)
Function for calculating scalar products given the information about the basis and the scalar product matrix in the basis. The Col amps are first decomposed using decompose, and then squared using the scalar product matrix P spm. An orthogonal scalar product matrix is assumed.

void scalar product matrix()
Calculates the scalar product matrix assuming the basis to be orthogonal. Calculates both the double (d spm) and the Polynomial (P spm) matrices and saves to default filenames. Function for calculating scalar products given the information about the basis and the scalar product matrix in numerical form. The Col amps are first decomposed using decompose, and then squared using diagonal d spm. For Orthogonal basis an orthogonal scalar product matrix is assumed. Calculates the scalar product between decomposed amplitudes v1, v2 using the diagonal d spm diagonal numerical scalar product matrix. The vectors needs to be expressed in the basis contained in cb, i.e., the decomposition has to be known.