# 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.

`global_def`

: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.`functions`

: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.

`collier`

: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`

.`phase_space`

: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.`openloops`

: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.`olinterface`

: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.`{pg}_mat_el`

: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.`{pg}`

: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

`mat_el`

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

`user`

: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.`vegas`

: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`

.`integrands`

: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.`mcmule`

:This is the main program, but actually does little else than read the inputs and call

`vegas`

with a function provided by`integrands`

.`test`

: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.

When started,

`mcmule`

reads options from`stdin`

as specified in Table 1.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.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\).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.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.

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`

.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.This triggers the computation of the matrix element and the assembly of the full weight.

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.

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

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.

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.