Structure of McMule

McMule is written in Fortran 95 with helper and analysis tools written in python [1]. The code is written with two kinds of applications in mind. First, several processes are implemented, some at NLO, some at NNLO. For these, the user can define an arbitrary (infrared safe), fully differential observable and compute cross sections and distributions. Second, the program is set up such that additional processes can be implemented by supplying the relevant matrix elements.

Modular structure of the code

McMule consists of several modules with a simple, mostly hierarchic structure. The relation between the most important Fortran modules is depicted in Figure 5. A solid arrow indicates “using” the full module, whereas a dashed arrow is indicative of partial use. In what follows we give a brief description of the various modules and mention some variables that play a prominent role in the interplay between the modules.

graph TD openloops(openloops) --> olinterface global_def --> olinterface collier(collier) --> functions global_def --> functions functions --> phase_space olinterface --> pgmatel subgraph " " pgmatel["{pg}_mat_el"] --> pg["{pg}"] end pg --> mat_el mat_el --> integrands functions --> user phase_space --> integrands user --> integrands integrands --> test integrands --> mcmule vegas --> test vegas --> mcmule phase_space -.->|"ksoft"| pgmatel user -.-> |"metadata"| vegas vegas -.-> |"bin_it"| integrands

Figure 5 The structure of McMule


This module simply provides some parameters such as fermion masses that are needed throughout the code. It also defines real(kind=prec) as a generic type for the precision used. [2] Currently, this simply corresponds to double precision.


This module is a library of basic functions that are needed at various points in the code. This includes dot products, eikonal factors, the integrated eikonal, and an interface for scalar integral functions among others.


This is an external module [3, 4, 5, 6]. It will be linked to McMule during compilation and provides the numerical evaluations of the scalar and in some cases tensor integral functions in functions.


The routines for generating phase-space points and their weights are collected in this module. Phase-space routines ending with FKS are prepared for the FKS subtraction procedure with a single unresolved photon. In the weight of such routines a factor \(\xi_1\) is omitted to allow the implementation of the distributions in the FKS method. This corresponds to a global variable xiout1. This factor has to be included in the integrand of the module integrands. Also the variable ksoft1 is provided that corresponds to the photon momentum without the (vanishing) energy factor \(\xi_1\). Routines ending with FKSS are routines with two unresolved photons. Correspondingly, a factor \(\xi_1\,\xi_2\) is missing in the weight and xiout1 and xiout2, as well as ksoft1 and ksoft2 are provided. To ensure numerical stability it is often required to tune the phase-space routine to a particular kinematic situation.


This is the external OpenLoops library [1, 2] that we use for some real-virtual matrix elements. It is pulled as a git submodule and linked to McMule during compilation.


This connects openloops to the rest of McMule by initialising OpenLoops for the process under consideration and converting to and from the OpenLoops conventions which are slightly different than the ones used by McMule.


Matrix elements are grouped into process groups such as muon decay (mudec) or \(\mu\)-\(e\) and \(\mu\)-\(p\) scattering (mue). Each process group contains a mat_el module that provides all matrix elements for its group. Simple matrix elements are coded directly in this module. More complicated results are imported from sub-modules not shown in Figure 5. A matrix element starting with P contains a polarised initial state. A matrix element ending in av is averaged over a neutrino pair in the final state.


In this module the soft limits of all applicable matrix elements of a process group are provided to allow for the soft subtractions required in the FKS scheme. These limits are simply the eikonal factor evaluated with ksoft from phase_space times the reduced matrix element, provided through mat_el.

This module also functions as the interface of the process group, exposing all necessary functions that are imported by


which collects all matrix elements as well as their particle labelling or PID.


For a user of the code who wants to run for an already implemented process, this is the only relevant module. At the beginning of the module, the user has to specify the number of quantities to be computed, nr_q, the number of bins in the histogram, nr_bins, as well as their lower and upper boundaries, min_val and max_val. The last three quantities are arrays of length nr_q. The quantities themselves, i.e. the measurement function, is to be defined by the user in terms of the momenta of the particles in quant(). Cuts can be applied by setting the logical variable pass_cut to false [3]. Some auxiliary functions like (pseudo)rapidity, transverse momentum etc. are predefined in functions. Each quantity has to be given a name through the array names.

Further, user contains a subroutine called inituser(). This allows the user to read additional input at runtime, for example which of multiple cuts should be calculated. It also allows the user to print some information on the configuration implemented. Needless to say that it is good idea to do this for documentation purposes.


As the name suggests this module contains the adaptive Monte Carlo routine vegas [15] . The binning routine bin_it is also in this module, hence the need for the binning metadata, i.e. the number of bins and histograms (nr_bins and nr_q, respectively) as well as their bounds (min_val and max_val) and names, from user.


In this module the functions that are to be integrated by vegas are coded. There are three types of integrands: non-subtracted, single-subtracted, and double-subtracted integrands, corresponding to the three parts of the \({\rm FKS}^2\) scheme [8, 25]. The matrix elements to be evaluated and the phase-space routines used are set using function pointers through a subroutine initpiece. The factors \(\xi_i\) that were omitted in the phase-space weight have to be included here for the single- and double-subtracted integrands.


This is the main program, but actually does little else than read the inputs and call vegas with a function provided by integrands.


For developing purposes, a separate main program exists that is used to validate the code after each change. Reference values for matrix elements and results of short integrations are stored here and compared against.

The library of matrix elements deserves a few comments. As matrix elements quickly become very large, we store them separately from the main code. This makes it also easy to extend the program by minimising the code that needs to be changed.

We group matrix elements into process groups, generic processes, and generic pieces as indicated in Appendix Available processes and which_piece. The generic process is a prototype for the physical process such as \(\ell p\to\ell p\) where the flavour of the lepton \(\ell\) is left open. The generic piece describes a part of the calculation such as the real or virtual corrections, i.e. the different pieces of (8) (or correspondingly (14) at NNLO), that themselves may be further subdivided as is convenient. In particular, in some cases a generic piece is split into various partitions (cf. Section Phase-space generation for details on why that is important).

What happens when running

In the following we discuss what happens behind the scene when asking McMule to perform the calculation of m2enng0 in Section Simple runs at LO.

  1. When started, mcmule reads options from stdin as specified in Table 1.

  2. Once McMule knows its configuration it associates the numerical values of the masses, as specified through flavour. In particular, we set the generic masses Mm and Me to Mtau and Mel. This is done in init_flavour(), defined in global_def. For other processes this might also involve setting e.g. centre-of-mass energies scms to default values.

  3. Next, the function to be integrated by vegas is determined. This is a function stored in integrands. There are basically three types of integrands: a standard, non-subtracted integrand, sigma_0, a single-subtracted integrand needed beyond LO, sigma_1, and a double-subtracted integrand needed beyond NLO, sigma_2. Which integrand is needed and what matrix elements and phase-space it depends on is determined by calling the function init_piece which uses the variable which_piece to point function pointers at the necessary procedures. For our LO case, init_piece sets the integrand to sigma_0 and fixes the dimension of the integration to \({\tt ndim}=8\).

  4. The matrix element pointer is assigned to the matrix element that needs to be called, Pm2enngAV(q1,n1,q2,q3,q4,q5). The name of the function suggests we compute \(\mu(q_1,n_1)\to [\nu(q3) \bar\nu(q4)]e(q_2) \gamma(q_5)\) with the polarisation vector n1 of the initial lepton. Even though we average over the neutrinos, their momenta are still given for completeness.

  5. The interplay between the function sigma_0(x,wgt,ndim) and vegas is as usual, through an array of random numbers x of length ndim that corresponds to the dimension of the integration. In addition there is the vegas weight of the event, wgt due to the Jacobian introduced by the importance sampling. The function sigma_0 simply evaluates the complete weight wg of a particular event by combining wgt with the matrix element supplemented by symmetry, flux, and phase-space factors.

    1. In a first step a phase-space routine of phase_space is called. For our LO calculation, init_piece pointed a pointer to the phase-space routine psd5_25(), a phase-space routine optimised for radiative lepton decays (cf. Section Phase-space generation). This will be called as a first step in the integrand to generate the momenta with correct masses as well as the phase-space weight weight.

    2. With these momenta the observables to be computed are evaluated with a call to quant(). If one of them passes the cuts, the variable cuts is set to true.

    3. This triggers the computation of the matrix element and the assembly of the full weight.

    4. In a last step, the routine bin_it, stored in vegas, is called to put the weight into the correct bins of the various distributions. If the variable under- or overshoots the bounds specified by min_val and max_val, the event is placed into dedicated, infinitely big under- and overflow bins.

    These steps are done for all events and those after pre-conditioning are used to obtain the final distributions.

  6. After preconditioning the state of the integrator is reset, as is usual.

  7. During the main run, the code generates a statefile from which the full state of the integrator can be reconstructed should the integration be interrupted (cf. Section Differential distributions and intermediary state files for details). This makes the statefile ideal to also store results in a compact format.

  8. The value and error estimate of the integration is printed to stdout.

To analyse these results, we provide a python tool pymule, additionally to the main code for McMule. pymule uses numpy [26] for data storage and matplotlib for plotting [13]. While pymule works with any python interpreter, IPython [22] is recommended. We will encountered pymule in Section Analysing the results when we discuss how to use it to analyse results. A full list of functions provided can be found in Appendix pymule user guide.